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NOMENCLATURE 


C  Capacitance 

g  Acceleration  due  to  gravity 

h  Film  coefficient 

k  Thermal  conductivity 

L  (1)  Distance  from  surface  to  center  of  node,  also 
(2)  thickness  of  insulated  slab 

T  Transform  variable 

t  Temperature  of  body  under  consideration 

tj.  Temperature  of  bounding  fluid 

X  Distance  from  surface  to  point  under  consideration  in 
analytical  solutions 

a  Thermal  diffusivity 

AX  Thickness  of  node 

p  Density 

T  Time 

Tj  Time  at  end  of  ramp 
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FOREWORD 


This  report  documents  an  investigation  to  fill  an  existing 
gap  between  the  theory  and  the  application  of  finite  difference 
methods  to  transient  thermal  analyses.  A  methodology  is  developed 
that  provides  a  means  of  predicting  and  modifying  the  analytical 
error  associated  with  thermal  response  problems. 

The  work  was  performed  at  the  University  of  Nevada  during  the 
period  from  September  1966  through  September  1968  under  two 
separate  contracts  with  the  Naval  Weapons  Center,  China  Lake, 
California.  Contract  N60530-67-C-0051  for  the  period  from 
September  1966  to  September  1967  was  funded  by  the  Bureau  of 
Naval  Weapons  WepTask  RMMO-42-008/21671/F009-09-01  under  the 
cognizance  of  W.K.  Baker.  Contract  N60530-67-C-1278  for  the 
period  from  September  1967  to  September  1968  was  funded  by  WepTask 
A32~320“008/216-l/F008-09-01  under  the  cognizance  of  W.C.  Volz. 

The  technical  administrator  of  both  contracts  for  the  Naval 
Weapons  Center  was  L.D.  Schultz, 

Technical  reviewers  for  this  report  were  Dr.  William  H. 
Thielbahr,  of  the  Propulsion  Development  Department,  and  Leo  D. 
Schultz,  of  the  Weapons  Development  Department.  Appendixes  B, 

C,  and  D  containing  bibliographic  information  are  published  as 
received  under  contract;  they  have  not  been  edited  at  this  Center. 
This  report  is  released  at  the  working  level  for  information  only. 
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INTRODUCTION 


Uncritical  application  of  finite-difference  procedures  for  solving 
transient  aerodynamic  heat-transfer  problems  has  required  10  to  15 
hours  on  the  computer  for  each  hour  of  actual  flight.  In  addition,  the 
complexity  of  problem,  e.g.,  the  variety  in  geometry,  materials  of 
construction,  initial  conditions,  and  boundary  conditions,  is  such  that 
an  estimate  of  the  error  inherent  in  the  analysis  is  difficult. 
Consequently,  the  introduction  of  errors  of  an  intolerable  magnitude 
into  the  problem’s  solution  by  adjustment  of  the  spatial  increment  and/ 
or  time  increment  to  minimise  the  excessive  costs  caused  by  long 
computer  runs  is  difficult  to  avoid  and  even  more  difficult  to  evaluate. 
Of  the  many  sources  of  potential  significant  error,  such  as  inaccurate 
flight  data,  atmospheric  conditions,  aerodynamic  heat  transfer 
coefficient,  and  thermal  properties,  the  error  considered  here  is  the 
analytical  error  due  to  the  nature  of  the  finite-difference  approxima¬ 
tion  of  the  partial  differential  equations  of  transient  thermal 
response.  Because  this  error  can  be  as  great  as  20%  of  the  surface 
temperature  rise,  the  necessity  of  maintaining  the  analytical  error 
within  acceptable  limits  is  apparent. 

The  purpose  of  the  investigation  is  to  fill  the  existing  gap 
between  theory  and  practice  by  developing  a  methodology  that  provides 
the  engineer  with  a  means  for  predicting  the  analytical  error 
associated  with  a  specific  thermal  response  problem.  This  ability 
allows  the  selection  of  spatial  and  time  increments  such  that  the 
minimum  computation  time  is  assured  for  a  predetermined  allowable 
analytical  error  limit.  The  benefits  are  two-fold:  First,  the 
analytical  error  can  be  controlled  within  limits  governed  by  the 
problem  under  consideration;  and  second,  the  solution  is  obtained  for 
the  least  cost. 

A  brief  review  of  the  implicit  and  explicit  forms  of  the  finite 
difference  method  for  handling  complex  thermal  transients  would  be  in 
order  at  this  point.  In  both  methods,  a  system  of  equations  is 
written  describing  the  heat  transfer  processes  taking  place  in  a 
geometry  which  has  been  divided  into  discrete  nodes.  For  each  node,  an 
energy  balance  is  written:  The  sum  of  all  forms  of  energy  crossing 
the  node  boundaries  are  equated  to  the  time  rate  of  change  of  the 
heat  capacity  of  the  node.  Thus,  for  both  approaches,  explicit  and 
implicit,  a  transient  thermal  problem  generates  a  large  number  of 
equations,  specifically,  one  equation  for  each  node  in  the  system. 
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The  set  of  these  equations  is  to  be  solved  at  every  time  step  for  the 
temperature  of  each  node  in  the  system. 

It  is  at  this  point  that  the  difference  between  implicit  and 
explicit  becomes  apparent.  In  the  explicit  or  forward-difference 
method,  the  nodal  temperature  and  the  time  rate  of  change  of  the  node 
capacitance  are  referenced  to  the  beginning  of  each  time  step.  Thus, 
knowing  all  of  the  temperatures  at  the  beginning  of  a  time  step  for 
the  entire  system  of  nodes,  one  can  predict  the  temperature  each  node 
will  reach  at  the  end  of  the  time  step.  In  this  manner,  each  equation 
has  a  single  unknown,  the  node  temperature  at  the  end  of  the  time  step. 
Thence  comes  the  nomenclature  explicit:  each  equation  can  be  solved 
explicitly  for  its  single  unknown  temperature. 

The  implicit  approach  expresses  the  nodal  temperatures  and 
references  the  time  rate  of  change  of  the  node  capacitance  to  the  end 
of  each  time  step.  This  approach  also  results  in  a  system  of  equations, 
one  for  each  node  in  the  system;  but  the  individual  equations  in  the 
system  may  contain  several  unknowns.  Thus,  the  system  is  no  longer 
explicit  but  is  an  implicit  system  in  which  the  entire  set  of  equations 
must  be  solved  simultaneously. 

In  brief,  the  limitation  for  the  explicit  method  is  in  the  length 
of  time  step  which  may  be  taken  before  instability  sets  in.  For 
stability,  the  length  of  the  time  step  is  a  function  of  the  thinnest 
dimension  of  any  of  the  nodes  in  the  system  and,  for  practical  problems, 
can  be  as  small  as  thousandths  of  a  second.  This  renuirsment  can  lead 
to  an  exceedingly  large  number  of  time  steps  to  solve  practical 
transient  problems  and  can  take  an  excessive  amount  of  computer  time. 

On  the  other  hand,  the  implicit  method  has  no  limiting  time  step 
and  has  maximum  stability  compared  to  any  of  the  other  methods.  A 
thorough  discussion  of  the  stability  characteristics  of  the  implicit 
and  explicit  methods  is  presented  in  the  paper  titled  The  Stability 
of  Three  Finite  Difference  Methods  of  Solving  for  Transient 
Temperatures  by  G.R.  Gaumer.  (See  Entry  41,  Appendix  C) . 

In  the  one-dimensional  problem  shown  by  Fig.  1,  the  surface 
undergoes  some  form  of  a  boundary  condition  change,  and  an  energy 
balance  is  written  for  each  of  the  nodes.  The  resulting  system  of 
equations  in  which  all  nodal  temperatures  are  referenced  to  the  end  of 
the  time  step  will  form  the  implicit  system.  The  next  step  is  to  solve 
the  resulting  matrix  for  the  nodal  temperatures  either  by  a  direct 
method  such  as  the  tridiagonal  matrix  solution  or  by  an  iterative 
method  such  as  Gauss-Seidel .  After  determining  the  nodal  temperatures 
at  the  end  of  a  time  step,  a  new  set  of  equations  is  generated  using 
the  same  principles  and  arriving  at  a  second  matrix  which,  in  turn, 
will  determine  the  temperatures  at  the  end  of  the  second  time  step. 

This  process  is  repeated  for  each  succeeding  time  step.  In  making 
comparisons  of  the  effectiveness  of  several  techniques  for  solving 
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the  matrix,  referred  to  previously,  it  was  necessary 
classical  analytical  solution  for  an  infinite  plate, 
given  in  Appendix  A  for  the  convenience  of  the  user, 
eigenvalue  generator. 
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FIG.  1.  Nodal  Network  for  an  Insulated  Plate. 

To  improve  accuracy  and  to  provide  rapid  solutions,  some  form  of 
compromise  is  necessary  for  intelligent  use  of  the  finite-difference 
method.  In  general  but  not  under  all  conditions,  the  finer  the  sub¬ 
division  of  time  steps  in  the  time  network  or  the  finer  the  sub¬ 
division  of  geometrical  spaces,  i.e.,  the  thickness  of  each  layer,  the 
closer  the  numerical  analysis  will  come  to  approximating  the  analytical 
solution.  However,  the  finer  time  network  and  the  finer  the  spatial 
network  for  any  given  problem,  the  longer  will  be  the  computation  time. 
In  the  explicit  method  round-off  errors  can  become  significant.  In  the 
implicit  method,  however,  round-off  errors  are  fairly  well  confined 
to  an  individual  time  step,  and  any  error  that  is  transmitted  to  the 
next  time  step  is  diffused  among  the  nodes.  If  the  user  knows  the 
effect  a  given  spatial  network  will  have  in  terms  of  spatial-truncation 
error  and  time-truncation  error,  a  much  more  effective  use  of  computer 
time  may  be  made.  Furthermore,  acceleration  techniques  may  be  used  to 
shorten  the  computer  time  needed.  Round-off  and  truncation  errors 
as  well  as  acceleration  techniaues  are  discussed  in  Appendix  B. 

So  that  one  can  gain  insight  into  the  nature  of  the  truncation 
errors,  Fig.  2  shows  a  true  temperature-distance  plot  for  two  adjacent 
nodes.  Tangent  to  the  true  temperature  curve  are  two  lines,  True 
Slope  1-2  and  True  Slope  2-3.  These  lines  represent  the  true  slope  of 
the  temperature  curve  at  the  midpoint  betweer  node  one  and  node  two  and 
at  the  midpoint  between  node  two  and  node  three.  Compare  these  true 
slopes  with  the  linear  approximations  that  are  shown  as  dashed  lines  in 
Figure  2;  one  sees  that  spatial-truncation  error  consists  of  the 
difference  between  the  true  slope  and  the  linear  approximation  of  the 
slope.  Using  this  linear  approximation  of  the  temperature  gradient 
results  in  an  error  for  the  rate  of  change  of  energy  stored  in  the 
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node.  This,  in  turn,  causes  an  error  in  the  prediction  of  the  nodal 
temperature  at  the  end  of  the  time  step. 


It  has  been  found  that  by  using  the  Fourier  number  and  the  Biot 
number,  a  correlation  of  "he  spatial  truncation  error  may  be  made 
(Entry  7,  Appendix  C) .  In  like  manner,  time  truncation  errors  may  also 
be  correlated. 

Commonlv,  the  error  analysis  of  transient  problems  is  based  on 
step-function  boundary  conditions;  i.e.,  the  driving  potential  undergoes 
an  instantaneous  change  at  the  beginning  of  the  transient.  In  actual 
aerodynamic  heat  transfer  problems,  step  changes  at  the  boundary  are 
seldom  realized;  changes  occur  over  time  periods  of  significant 
duration.  For  this  reason,  ramp  functions  approximate  the  true 
boundary  condition  changes  of  practical  aerodvnamic  problems  more 
realistically  than  step  functions  do.  Therefore,  after  the  basic 
analytical  error  analysis  is  developed  for  step  function  changes  at  the 
boundary,  the  analysis  is  extended  to  determine  the  effect  of  duration 
and  slope  of  ramp  function  boundary  conditions  on  the  spatial  truncation 
error , 

The  summarv  of  a  literature  survey,  made  in  the  initial  months  of 
this  research  effort,  is  included  in  this  report  as  Appendix  B.  This 
literature  survey  compares  a  number  of  papers  concerned  with  the 
various  techniques  for  establishing  iterative  solutions,  techniques  for 
determining  the  truncation  errors,  both  spatial  and  time,  and  techniques 
for  accelerating  iterative  solutions.  Appendix  C  is  a  bibliography  on 
the  general  subject  of  numerical  solution  for  transient  heat  transfer 
problems;  Appendix  D  contains  the  same  bibliography  organized  by 
categories.  Aprendix  E  is  a  derivation  of  the  Laplace  transform 
solution  of  the  one-dimensional  Fourier  conduction  equation  for  a 
semi-infinite  solid  (and  for  an  insulated  slab)  with  a  ramp-function 
boundary  condition. 
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MATRIX  SOLUTION  METHODS 


Because  the  transient  heat-transfer  problems  in  this  study  are 
one-dimensional,  a  variety  of  solutions  to  the  resulting  matrix  are 
available.  For  example,  an  iterative  methods— -Gauss-Seidel — and  the 
various  acceleration  routines  which  may  be  used  with  it;  the  Runge- 
Kutta  method;  and  a  direct  method  utilizing  the  tridiagonal  form  of  a 
matrix  (TRIDAG)  are  available.  The  direct  method  bypasses  the 
difficulty  of  convergence  errors  found  in  iterative  solutions.  As  an 
example  J.O.  Wilkes  of  the  University  of  Michigan  reports  on  a 
solidification  problem  involving  11  nodes  and  17  time  steps  for  which 
the  IBM  7090  has  an  execution  time  of  1.8  seconds.  (See  Entry  110, 
Appendix  C.)  Comparable  problems  using  an  iterative  solution  took  an 
average  of  18.0  seconds  for  execution. 

For  nonlinear  problems,  a  simple  direct  solution  such  as  TRIDAG 
is  not  applicable;  therefore,  a  more  general  approach  such  as  the 
accelerated  Gauss-Seidel  using  a  constant  acceleration  factor  (usually 
referred  to  as  over-relaxation)  may  be  used.  An  alternative  is  to  use 
Gauss-Seidel  with  an  acceleration  technique  such  as  the  Wegstein  or 
the  Steffensen  methods. 

A  study  of  Steffensen's  accelerating  technique  applied  to  a  simple 
algorithm  indicated  possibilities.  However,  when  this  technique  was 
used  to  accelerate  the  solution  of  a  20-node  problem,  the  results  were 
either  the  same  as  could  be  expected  from  the  regular  Gauss-Seidel  or 
greater.  The  system  did  not  seem  to  settle  down  even  though  the 
frequency  of  application  of  the  technique  was  varied. 

The  Wegstein  method  requires  more  computer  storage  space  and  more 
computer  operations  than  does  the  Steffensen  method,  but  the  Wegstein 
method  effectively  reduces  the  number  of  iterations  needed  for 
convergence.  Appendix  A  presents  Wegstein' s  basic  equations  and 
incorporates  the  Wegstein  method  into  the  ONE-D  program  for  use  with  the 
Gauss-Seidel  iteration  technique.  A  FORTRAN  IV  listing  of  this  program 
is  given  in  Appendix  A. 

A  note  of  caution:  some  acceleration  techniques  such  as  Steffensen 
or  Wegstein  work  very  well  on  the  algorithm  type  of  iterative  solution, 
but  when  used  on  a  system  of  equations,  these  techniques  can  actually 
slow  down  the  rate  of  convergence.  For  comparison  purposes,  a  20-node 
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one-dimensional  problem  having  a  step  input  in  the  bounding  fluid 
temperature  was  devised.  Straight  Gauss-Seidel  using  a  1/100-degree 
convergence  limit  took  46  iterations  to  determine  the  nodal  temperatures 
at  the  end  of  the  first  time  step.  An  unrestrained  acceleration 
technique  such  as  Wegstein's  applied  every  third  iteration  caused  such 
divergence  that  the  automatic  program  stop  of  100  iterations  was 
reached.  This  difficulty  was  corrected  by  restricting  the  Wegstein 
technique  only  to  positive  values  of  acceleration.  Only  17  Iterations 
were  required  for  convergence  in  a  subsequent  computer  run  using  a 
combination  of  10  initial  iterations  of  Gauss-Seidel  without  over¬ 
relaxation,  followed  by  the  Wegstein  acceleration  method  limited  to 
positive  values  only,  and  then  three  additional  iterations  of  standard 
Gauss-Seidel.  As  a  further  experiment,  a  run  with  10  Gauss-Seidel 
over-relaxed  iterations,  followed  by  a  Wegstein  acceleration  limited 
to  positive  values  only,  followed  by  standard  Gauss-Seidel  alternating 
with  the  Wegstein  method  took  20  iterations  to  converge  the  same 
problem.  Mixing  combinations  of  various  acceleration  techniques  should 
be  approached  with  cautior  because  some  combinations  actually  increase 
the  number  of  iterations  over  that  required  by  a  single  technique. 


SPATIAL-,  TIME-,  AND  TOTAL-TRUNCATION  ERRORS 


One  of  the  most  confusing  problems  facing  the  engineer  in  attempt¬ 
ing  to  utilize  finite  differences  to  solve  transient  heat-transfer 
problems  is  the  determination  of  the  spatial-network  grid  and  the  time 
interval  to  take  in  establishing  a  time  network  for  dealing  with  a 
transient.  In  an  effort  to  quantize  this  particular  problem,  the 
available  background  material  was  gleaned  from  the  literature.  (See 
Entry  7,  9  and  44  in  Appendix  B.)  Truncation  errors  are  caused  by  the 
strong  second  derivatives  inherent  in  the  beginning  of  the  transient 
and  near  the  surface  of  the  geometrv.  The  curves  that  are  available 
have  been  found  to  be  quite  successful  in  predicting  spatial-  and 
time-truncation  errors  in  problems  consisting  of  homogeneous  materials 
having  uniform  thicknesses  of  the  individual  slices.  Some  of  these 
curves  are  given  in  Fig.  3,  4  and  5  in  which  the  abscissa  for  the 
truncation-error  curves  are  the  Fourier  number  (the  thermal  diffusivitv 
times  the  time  span  from  the  beginning  of  the  transient  up  to  the 
particular  instant  for  which  an  error  evaluation  is  desired,  divided  by 
the  square  of  the  distance  from  the  surface  to  where  the  error  is  being 
evaluated).  This  family  of  curves  is  correlated  by  the  Biot  number  (the 
film  coefficient  at  the  surface  of  the  geometry,  times  the  distance  from 
the  surface  down  to  the  point  in  question,  divided  by  the  thermal 
conductivity  of  the  material) .  The  measurement  of  the  error  itself  is 
done  in  terms  of  the  step-function  temperature  rise  in  the  bounding 
fluid  at  the  beginning  of  the  transient.  It  is  of  interest  to  note  in 
Fig.  3  that  when  the  thickness  of  the  slice  is  reduced,  which  is 
proportional  to  L  for  the  first  node,  the  Fourier  number  and  the  Biot 
number  are  both  affected  such  that  the  spatial-truncation  error  is 
reduced.  However,  it  may  be  noted  also  that  it  is  possible  in  a 
Fourier  number  range  from  approximately  0  to  1.0  to  reduce  the  thickness 
of  a  slice  and  have  the  spatial-truncation  error  increase. 

With  respect  to  variations  in  Biot  number,  one  can  easily  see 
that  the  worst  possible  case  occurs  when  the  film  coefficient  at  the 
surface  is  infinite,  and  the  errors  can  range  as  high  as  13%  of  the 
temperature  rise  in  the  bounding  fluid.  However,  with  decreasing  values 
of  the  film  coefficient,  reflected  in  decreasing  values  of  Biot  number, 
the  spatial-truncation  errors  are  in  turn  decreased.  As  a  result  it  is 
possible  for  certain  classes  of  problems,  that  the  spatial-truncation 
error  will  never  get  above  1  or  2%  regardless  of  the  thickness  of  the 
slice. 
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FIG.  3.  Spatial-Truncation  Error  for  Nodes  Adjacent  to  the  Surface. 

In  like  manner,  time-truncation  errors  can  be  evaluated  using  the 
curves  shown  in  Fig.  4.  Once  again,  the  Fourier  number  is  used  as  the 
abcissa  for  the  curves,  the  Biot  numbers  for  the  family  of  curves,  and 
the  errors  are  expressed  on  the  ordinate  as  a  percentage  of  the  step- 
function  temperature  change  in  the  bounding  fluid  at  the  beginning  of 
the  transient.  Once  again,  we  see  the  effect  of  the  Biot  number;  i.e., 
as  the  Biot  number  becomes  smaller,  the  errors  become  smaller.  In  a 
final  evaluation,  one  sees  in  Fig.  5  that  the  errors  are  accumulative 
and  that  the  total  error  under  the  worst  possible  set  of  circumstances 
could  be  as  high  as  20%  of  the  step  function.  One  important  difference 
between  the  time-truncation-error  curves  and  the  spatial-truncation- 
error  curves  should  be  realized.  The  time  used  in  the  Fourier  number 
for  the  time-truncation  error  curves  is  the  time  interval  between  the 
beginning  of  the  transient  and  the  first  evaluation  of  temperature  in 
the  problem;  i.e.,  the  first  time  step.  The  time  in  the  Fourier  number 
for  the  spatial-truncation  error  curves  refers  to  the  total  time 
interval  from  the  beginning  of  the  transient  to  the  time  at  which  an 
error  evaluation  is  made,  regardless  of  how  many  time  steps  elapsed. 
This  means  that  the  time-  and  spatial-truncation  errors  may  be  super¬ 
imposed  only  for  the  first  time  step.  It  is  also  important  to  realize 
that  the  distance,  L,  in  the  dimensionless  groups  of  Fourier  and  Biot 
numbers  would  be  half  the  thickness  of  the  node  slice  when  considering 
the  first  node  adjacent  to  the  surface  and,  subsequently,  would  be 
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FIG.  4.  Time-Truncation  Error  Versus  Fourier  Number  for  the 
First  Node  in  Homogeneous  Slab. 
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FIG.  5.  Truncation  Error  Versus  Fourier  Number  for  First 
Node- First  Time  Step  in  Homogeneous  Slab. 
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three  halves  the  thickness  of  one  node  when  considering  the  errors  for 
the  second  node  from  the  surface,  five  halves  of  the  node  thickness 
when  dealing  with  the  third  node  from  the  surface,  etc. 

When  dealing  with  problems  in  which  layers  of  nonhomogeneous 
materials  are  being  considered,  or  in  which  a  ramp  function  instead 
of  a  step  function  occurs  on  the  surface,  or  when  a  nonflat  plate 
geometry  is  encountered,  the  truncation  errors  are  no  longer  exact. 
However,  in  the  practical  sense,  valuable  guidance  may  be  obtained 
from  these  curves. 

In  considering  the  effect  of  successive  nodes  on  spatial-truncation 
errors  going  into  depth  from  the  surface  in  a  homogeneous  flat  plate, 
the  trend  is,  in  all  cases,  for  the  greatest  truncation  error  to  occur 
in  the  node  adjacent  to  the  surface;  and  with  each  succeeding  node, 
the  truncation  error  is  reduced  (Fig.  6).  Generally,  this  reduction 
in  spatial-truncation  error  is  quite  drastic;  by  the  time  the  third 
or  fourth  node  is  reached,  the  truncation  errors  are  of  the  order  of 
magnitude  of  1%  or  less.  One  can  see  that  time-truncation  errors  are 
rapidly  reduced  if  two  or  more  time  steps  are  utilized  to  complete  the 
first  time  interval  (Fig.  7,  8,  and  9).  Also  the  time-truncation  error 
becomes  less  for  successive  nodes  in  depth  (Fig.  10) . 

It  is  concluded  that  if  the  first  time  step  is  sized  to  keep  the 
time-truncation  error  within  reasonable  limits  in  a  homogeneous 
geometry,  the  time-truncation  error  for  the  nodes  below  the  first 
node  will  automatically  have  lower  error  levels. 


FIG.  6.  Spatial-Truncation  Error  for  Second  Node  From 
Surface. 
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FIG.  7.  Time-Truncation  Error  for  Nodes  Adjacent  to  the 
Surface;  Two  Time-Step  Numerical  Solution. 
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FIG.  8.  Time-Truncation  Error  for  Nodes  Adjacent  to  the 
Surface;  Four  Time-Step  Numerical  Solution. 
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FIG.  9.  Time-Truncation  Error  for  Nodes  Adjacent  to  the 
Surface;  Ten  Time-Step  Numerical  Solution. 
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FIG.  10.  Time-Truncation  Error  for  Several  Nodes. 


NUMERICAL  ERRORS  IN  THIN-THICK  GEOMETRIES 

A  recurring  type  of  problem  of  considerable  interest  is  the 
transient  thermal  analysis  of  a  two-layer  geometry.  The  two  layers 
consist  of  a  thin,  highlv  conductive  layer,  such  as  aluminum,  exposed 
on  one  side  to  a  bounding  fluid  and  on  the  other  side  to  a  “.hick  layer 
of  low-conducting  material  acting  as  an  insulator.  The  questions  to 
be  resolved  are  What  are  the  numerical  errors  inherent  in  such  a 
geometrical  system?  and  What  procedure  might  be  followed  to  keep  the 
errors  to  an  acceptable  level? 

The  truncation-error  curves  were  prepared  in  much  the  same  way 
as  those  presented  by  Graybeal  (Entry  44,  Appendix  C) .  However,  a 
somewhat  different  scheme  was  employed  to  obtain  a  reference  solution. 
Whereas  an  analytical  solution  was  employed  directly  by  Graybeal  as 
an  errorless  reference,  such  a  solution  was  used  indirectly  here.  This 
indirect  reference  solution  was  taken  from  a  tabulation  of  an 
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analytical  solution  to  the  one— dimensional  conduction  problem  involving 
a  two-layer  plate. ^  The  tables  were  not  used  directly  because  the 
results  listed  for  each  set  of  parameters  present  a  limited  coverage 
of  the  time-space  field,  necessitated  by  a  desire  to  offer  solutions  for 
a  large  number  of  combinations  of  parameters,  yet  keep  the  report  from 
being  excessively  bulky.  To  obtain  a  usable  reference  solution,  a  set 
of  parameters  was  chosen  corresponding  to  one  section  of  the  table 
(a  =  0.05) ;  a  finite-difference  solution  was  prepared  using  these 
parameters  and  also  using  very  small  AT  and  AX.  This  solution  was  then 
compared  to  the  tabulated  solution  at  points  of  correspondence.  A 
truncation  error  of  less  than  0.12  (based  on  the  size  of  the  step 
change  of  temperature  in  the  bounding  fluid)  was  found  at  each  point 
checked,  and  this  finite-difference  solution  was  then  used  as  the 
standard  for  evaluating  truncation  error  in  other  finite-difference 
approximations  employing  larger  AT  and  AX.  First,  spatial-truncation 
error  only  was  introduced  by  increasing  AX  and  holding  AT  at  the 
original  small  value,  then  evaluating  spatial- truncation  error  for  a 
given  AX  by  comparison  with  the  reference  solution.  Then,  with  AX  held 
at  the  same  value  used  in  the  reference  solution,  AT  was  varied  to 
give  a  solution  containing  time-truncation  error  alone. 

The  problem  considered  here  is  one-dimensional  conduction  in  a  two- 
layer  composite  slab  heated  on  one  face  by  convection  from  a  fluid  wh^ch 
undergoes  a  step  change  of  temperature.  The  fluid  and  slab  are 
initially  at  the  same  uniform  temperature  ,  and  the  face  of  the  slab  not 
in  contact  with  the  fluid  is  perfectly  insulated.  All  physical  proper¬ 
ties  are  assumed  constant. 

Two  cases  were  studied:  in  one  ,  the  thickness  of  the  thin  layer 
next  to  the  fluid  was  5%  of  the  total  thickness  of  the  slab;  in  the 
other,  the  thickness  of  the  thin  layer  was  20%.  The  physical 
properties  of  this  layer  correspond  roughly  to  those  of  aluminum,  and 
the  properties  of  the  thick  layer  are  similar  to  those  of  some  common 
insulating  materials; 


Properties 

Thin  layer 

Thick  laver 

k,  BTU/ (hr-f t-°F) . 

.  100 

0.1 

p,  lb/(ft3) . 

.  200 

20 

C,  BTU/ (lb-°F) . 

.  0.25 

0.25 

'''Naval  Ordnance  Test  Station.  Temp&ratuve  Tables.  Part  7 ,  Vol.  1 
and  2.  Two-Layer  Plate ,  One-Space  Variable,  Linear,  by  H.  N.  Browne,  Jr. 
and  C.  J.  Thorn.  China  Lake,  Calif.  ,  NOTS,  1  March  1960.  (NAV0RD 
Report  5562,  Part  7;  NOTS  TP  2182,  Vol.  1  and  2.) 
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RESULTS 

The  spatial-truncation  errors  for  the  first  node  of  the  first  layer 
and  the  first  and  second  node  of  the  second  layer  for  the  thin-thick 
geometry  are  plotted  in  Fig.  11. 

FOURIER  NUMBER.  a2T/L2 

0.01  0.05  0.10  0.50  1.00  6.00 


FIG.  11.  Spatial-Truncation  Error  Vs.  Fourier  Number  for 
Selected  Nodes  in  Thin-Thick  Configurations. 


These  curves  contain  spatial-truncation  errors  for  the  worst  case 
considered:  where  the  thin  layer  is  represented  by  a  single  node,  and 
all  nodes  are  the  same  size.  The  maximum  spatial-truncation  error  in 
the  thin  layer  is  0.3%,  while  the  maximum  at  the  first  node  of  the 
second  layer  is  -2.3%.  This  drops  to  a  maximum  of  -1.0%  at  the  second 
node  of  the  second  layer. 

It  may  be  noted  that  the  errors  plotted  in  Fig.  11  are  much  smaller 
than  the  errors  encountered  in  the  curves  describing  truncation  errors 
in  homogeneous  slabs:  for  example,  a  film  coefficient  of  1,000  and  a 
thermal  conductivity  of  25. 

In  the  thin-thick  case,  approximations  of  real  problems  always 
results  in  very  small  Biot  numbers.  For  example,  a  1/8-inch  thick 
layer  of  stee]  could,  in  the  limit,  form  a  single  node  and  under 
extreme  circumstances,  give  rise  to  a  Biot  number  of  0.2  (h  -  1,000, 
k  -  25,  L  =  0.0625).  The  spatial-truncation  error  associated  with 
these  small  values  of  Biot  number  is  less  than  1%  (Fig.  3). 
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To  sum  up  then,  the  curve  in  Fig.  li  labeled  first  node  of  the 
first  layer  having  only  a  fraction  of  a  percent  error  certainly  fits  the 
trend  discovered  in  homogeneous  slabs.  In  a  similar  fashion  the  curve 
labeled  first  node  second  layer,  the  Biot  number  has  a  value  of  7.5 
based  upon  again  a  film  coefficient  of  50  and  the  thermal  conductivity 
of  0.1,  which  represents  the  insulation  involved.  Once  again,  comparing 
this  with  the  results  for  homogeneous  slabs  in  Fig.  3,  one  sees  that 
the  homogeneous  slab  would  have  spatial-truncation  errors  of  approxi¬ 
mately  10%  for  this  Biot  number;  whereas  in  the  thin-thick  case,  the 
maximum  error  was  under  2  1/2%.  This  is  a  logical  extension  of  the 
homogeneous  material  problem  because  the  aluminum  layer  on  the  surface 
acts  as  a  buffer  and  tends  to  protect  the  insulating  layer  from  the 
extreme  spatial-truncation  errors  that  would  occur  if  the  insulation 
were  directly  exposed  to  the  bounding  fluid.  With  reference  to  the 
third  curve  in  Fig.  11  (that  is,  the  second  node  of  the  second  layer 
curve),  its  relationship  to  the  f4rst  node  of  the  second  layer  is 
approximately  the  same  as  second  nodes  normally  have  in  homogeneous 
materials,  as  shown  in  Fig.  6.  Thus,  the  data  in  Fig.  11  follows 
the  same  trends  as  the  data  for  homogeneous  slabs. 

The  time-truncation  errors  for  the  thin-thick  class  of  problem  are 
plotted  in  Fig.  12.  The  ordinate  in  Fig.  12  is  the  usual  error 
measured  in  percent  of  the  step-function  temperature  change  in  the 
bounding  fluid,  but  the  abscissa  is  the  percent  of  the  boundary 
temperature  step-function  change  that  equals  the  temperature  change  in 
the  first  node  of  either  the  first  or  second  layer  that  occurs  in  the 
first  time  step.  The  time-truncation  errors  for  a  similar  case 
involving  a  homogeneous  slab  are  also  displayed  in  Fig.  12. 

It  is  surprising  to  find  that  the  thin-thick  layer  has  somewhat 
larger  time-truncation  errors  than  the  homogeneous  slab  and  at  the 
present  time,  no  explanation  is  available.  The  curves  in  Fig.  12, 
however,  provide  a  very  practical  guide  for  those  who  are  concerned 
with  time-truncation  errors  because  the  temperature  response  of  the 
first  node  of  either  the  thin  or  the  thick  layer  can  be  expressed  as  a 
percent  of  the  boundary  temperature  step-function  change.  A  reasonable 
criterion  for  limiting  the  time-function  error  at  the  beginning  of  a 
transient  is  to  ensure  that  the  temperature  change  of  the  first  node  is 
not  greater  than  20%  of  the  step-function  boundary  temperature  change 
during  the  first  time  step.  At  the  20%  level  or  below,  all  of  the  time- 
truncation  errors  are  less  than  2%. 
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FIG.  12.  Time  Truncation  Error  Vs.  Temperature  Rise  in  Designated 
Nodes  Expressed  as  a  Percent  of  Boundary  Temperature  Step  Function. 
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RAMP-FUNCTION  SPATIAL-TRUNCATION  ERRORS 


Thus  far,  only  errors  for  step-funccion  boundary  temperature 
changes  have  been  considered.  Since  the  ramp  function  is  a  closer 
approximation  to  the  actual  boundlng-fluid  temperature  variations,  such 
mathematical  models  of  transient  heat-transfer  situations  are  more 
realistic,  and  the  system  response  should  result  in  a  lower  spatial- 
truncation  error. 

As  a  starting  point,  two  approaches  were  used  to  determine  the 
true  temperature  for  use  as  a  comparison  basis  for  determining  spatial- 
truncation  errors.  The  first  approach  employs  successively  smaller 
and  smaller  time  steps  with  smaller  and  smaller  geometrical  divisions 
in  numerical  solution  for  a  transient  one-dimensional  homogeneous  heat 
transfer  problem.  Thus  an  asymptotic  approach  to  the  true  temperature 
distribution  was  obtained.  The  other  approach  uses  the  Laplace 
transform  method  in  a  classical  analytical  solution  of  the  one¬ 
dimensional  Fourier  transient-conduction  equation.  Both  approaches  are 
successful,  and  a  comparison  of  the  results  of  the  successive  approxi¬ 
mation  method  with  the  results  from  the  analytical  method  is  given  in 
Table  1. 


TABLE  1.  Comparisons  of  Numerical  &  Analytical  Solutions 
for  Ramp  Functions  Boundary  Conditions. 


These  data  were  generated  by  a  problem  having  the  following 
specifications:  Semi-infinite  solid,  2500°F/sec  on  bounding 
fluid.  AX  =  0.006  in;  AT  =  0.0002  sec;  h  =  10^;  a  =  0.2. 


Temperature,  °F 

Time, 

Sec 

Node  1 

Node  2 

Node  3 

Analytic 

1 

Numerical 

Analvtic 

Numerical 

Analytic 

Numerical 

0.02 

37.919 

36.787 

20.814 

20.178 

10.658 

10.391 

0.04 

82.373 

81.159 

54.649 

53.798 

35.057 

34.523 

0.06 

128.111 

126.861 

92.090 

91.140 

64.752 

64.076 

0.08 

174.516 

173.245 

131.449 

130.439 

97.400 

96.634 

0.10 

221.348 

220.348 

172,049 

170.997 

131.999 

131.169 
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The  Laplace  transform  approach  is  summarized  below;  the  complete 
derivation  is  given  in  Appendix  E. 

The  partial  differential  eouation  for  a  one-dimensional  transient 
problem  is 


3t 

3T 


(1) 


where 


t  =  temperature 


I 


! 

s 

i 


T  =  time 

The  initial  conditions  are  at  T  =  0,  t  =  0.  From  the  Laplace  transform 
o  f  Ecj .  1 , 


2, 

dX 


(2) 


A  general  solution  of  this  second  order  differential  equation  is 


5. 

r 

I 

1 


I 


T  =  A  exp ( \fs/a  X)  +  B  exp(-v§7ct  X)  (3) 

where  A  and  B  are  constants.  Since  a  semi-infinite  geometry  is  under 
consideration,  A  =  0.  Constant  B  mav  be  found  from  the  boundary 
condition. 


S 

The  bounding  fluid  temperature  change  is  a  ramp  function  (Fig.  13). 


FIC.  13.  Bounding  Fluid  Temperature  Change. 
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At  the  surface,  X  =»  0. 

-k  =  h[f(T)  -  t(0,  T) 3 


(4) 


dT 


If  the  Laplace  transform  of  Eq.  4  is  taken  and  if  —  (0,s)  and 
T(0,s)  are  evaluated  from  Eq.  3  the  constant  B  may  he  expressed  as 

ht . 


B  - 


kT. 


1  -  exp(-T^s) 

s2(x/57<*  +  h/k) 


Thus, 


ht^  | exp(-^/s/a  X)  jl  -  exp(-T^s)J 
1 


T  *  W7 


s  (\Zs/a  +  h/k) 
The  inverse  transform  of  Eq.  6  is 


t=-rr((T1+|-+  — ^  + 
Tx  M  1  2a  ah 


°  (s)  )  ertc(2  jrr) 


1  /  k\2  [h  '■  '2 

a  h 


t  = 


-C(x+^)  -'’[-(i^)2]!'  °  <  T<-h 


“(t)  T]  erfc(v§r + 1  y^) 


2a 


Si  +i 

ah  a  \h/ 


erfc 


X 

—  ■  7— w mctssm 

2Va  (T-Tj. 


5  (S' 


(5) 


(6) 


(7) 
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Note  chat  since  there  is  a  discontinuity  at  time  Ti,  there  are  two 
solutions:  one  for  T  greater  than  zero,  but  equal  to  or  less  than  T].; 
and  another  for  T  greater  than  Tj.  It  may  be  further  noted  that  at 
time  greater  than  Tj,  the  first  portion  of  the  equation  is  the  same  as 
the  equation  for  time  less  than  T^;  and  that  the  second  portion  of  the 
equation  is  similar  to  the  first  portion  except  T  is  replaced  by 
(T-  ?p) .  This,  in  effect,  applies  a  negative  slope  at  time  Tj  as  is 
illustr-ted  in  Fig.  14. 


TIME 


FIG.  14.  Synthesis  of  Analytical  Solution 


Figure  14  illustrates  the  concept  of  subtracting  At,  which  develops 
after  T^,  from  fluid  temperature  which  would  have  developed  after  T]_  if 
the  ramp  had  continued. 
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The  effect  is  that  the  first  portion  of  Eq.  8  will  produce  node 
temperatures  analogous  to  an  extension  of  the  ramp  after  Ti  and  the 
last  half  of  Eq.  8  will  produce  a  At  generated  in  the  nodes  after  Ti> 

In  similar  manner,  the  analytical  solution  for  the  case  of  a  ramp- 
function  temperature  change  in  the  bounding  fluid  adjacent  to  the 
surface  of  an  insulated  slab  was  derived.  The  final  equation  for  the 
insulated  slab  is 


A  complete  derivation  of  Eq.  9  is  given  in  Appendix  E. 

Computer  programs  were  coded  in  FORTRAN  for  both  of  the  ramp- 
function  analytical  solutions,  the  semi-infinite  solid  and  the  insulated 
slab.  FORTRAN  IV  listings  for  these  two  programs  are  included  in 
Appendix  A  of  this  report.  These  analytical  solutions  were  used  to 
calculate  the  temperatures  for  comparison  with  the  temperatures  for 
equivalent  ramp-function  problems  solved  by  numerical  methods.  The 
differences  between  results  of  these  two  methods  in  terms  of  temperature 
were  the  spatial-truncation  errors  sought.  The  FORTRAN  IV  listing  for 
the  numerical  method  computer  program — coded  such  that  either  step- 
function  temperature  rise  in  the  bounding  fluid  or  a  ramp-function 
temperature  rise  in  the  he  bounding  fluid  can  be  specified — is  given  in 
Appendix  A.  The  attempt  to  make  the  same  type  of  correlation  that  was 
successful  in  the  step-function-error  curve  analysis  failed,  as  shown  in 
Fig.  15.  That  is,  the  assumptions  that  the  errors  could  be  correlated 
by  the  overall  temperature  rise  in  the  fluid  and  that  a  further  correla¬ 
tion  would  be  available  based  upon  the  value  of  the  Biot  number  is 
incorrect;  the  Biot  number  does  not  correlate  the  spatial -truncation 
error  for  ramp-function  boundary  conditions.  Figure  16  shows  that  to 
correlate  various  ramp  functions  for  the  same  Biot  number  using  the 
overall  temperature  rise  in  the  bounding  fluid  as  the  basis  for  percent 
error  is  not  successful. 

A  correlation  using  the  current  fluid  temperature  rise  as  the  basis 
for  establishing  the  percent  spatial-truncation  error  causes  all  the 
ramp-function  errors  to  fall  on  the  same  curve  (Fig.  17).  The  data 
presented  in  Fig.  17  were  obtained  from  runs  using  a  long  ramp.  A 
long  ramp  is  one  in  which  the  maximum  spatial-truncation  error  occurs 
well  before  the  end  of  the  ramp  is  reached.  Of  interest  in  Fig.  17 
is  the  maximum  error  of  about  -9%  which  occurs  at  a  Fourier  number  of 
about  2.0.  These  values  are  to  be  compared  with  the  -13%  error  at  a 
Fourier  number  of  about  0.9  for  the  step  function  boundarv  condition. 
Maximum  error  is  reduced  because  of  the  reduced  rate  of  change  at  the 
model  boundary,  and  the  Fourier  number  is  increased  due  to  the  increased 
time  to  reach  terminal  fluid  temperatures. 
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Using  the  current  fluid  temperature  approach,  an  investigation  was 
made  in  which  the  ramp  temperature  rise  in  the  bounding  fluid  was 
completed  before  any  spatial-truncation  errors  could  develop  in  the 
first  node.  This  situation  is  termed  an  ultra-short  ramp.  The  ultra- 
short  ramp  curve  is  shown  in  Fig.  18;  the  dotted  line  represents  the 
equivalent  spatial-truncation  error  curve  for  a  step-function 
temperature  rise  in  the  bounding  fluid.  Again,  the  ramp-function 
error  curve  is  slightly  delayed  compared  to  the  step-function  error 
curve . 

A  further  study,  graphed  in  Fig.  19,  shows  the  result  when  the  end 
of  the  ramp  occurs  during  the  development  of  the  spatial-truncation 
error,  termed  a  short  ramp.  The  arrow  in  Fig.  19  indicates  the  time  of 
which  the  end  of  the  ramp  occurs;  after  the  end  of  the  ramp,  a  typical 
step-function  curve  develops.  Thus  Fig.  19  is  an  excellent  example  of 
how  the  nodal  system  shifts  from  a  ramp-function  error  response  to  a 
step-function  error  response  and  how  the  Fourier  number  for  the  maximum 
error  is  increased  due  to  the  ramp-function  effect. 

In  Fig.  20,  the  arrows  again  indicate  the  Fourier  number  corre¬ 
sponding  to  the  end  of  the  ramp  temperature  rise  in  the  bounding  fluid. 
For  these  curves  the  end  of  the  ramp  occurs  shortly  after  the  maximum 
spatial  truncation  error  occurs.  It  may  be  seen  that  after  the  end  of 
the  ramp,  the  individual  runs  exhibit  small  increases  in  error  showing 
a  residual  tendency  to  respond  in  a  step  function  fashion. 

These  combined  data  on  one  set  of  coordinates  are  shown  in  Fig.  21, 
the  dotted  line  represents  the  equivalent  step-function  error  curve  and 
the  solid  lines  represent  short-vamp  functions  and  long-vamp  functions. 
Of  interest  are  the  bumps  or  departures  (dashed  lines)  visible  in  the 
short-ramp  functions.  These  departures  represent  a  progression  of 
error  increases  that  occur  when  the  ramp  portion  of  the  rise  in  the 
bounding  temperature  ceases  before  the  maximum  spatial-truncation  error 
occurs.  At  the  end  of  the  ramp,  the  error  rapidly  climbs  toward  the 
levels  exhibited  by  a  step-function  rise  in  bounding  fluid  temperature. 
In  this  figure  the  long-ramp  function  (a  ramp  in  which  the  maximum 
spatial-truncation  error  occurs  while  the  ramp  is  still  in  progress) 
properly  defines  the  error  expressed  as  percent  of  current  temperature 
rise  in  the  bounding  fluid.  For  low  Biot  numbers,  it  makes  little 
practical  difference  whether  or  not  a  ramp  function  or  a  step  function 
occurs  in  the  bounding  fluid  because  spatial-truncation  errors  for  low 
Biot  numbers  are  small  in  either  case.  As  an  example,  the  maximum 
possible  spatial-truncation  error  is  2.2%  for  a  Biot  number  of  0.5. 

Whether  or  not  a  given  problem  has  a  short-  or  a  lo^g-ramp  bounding 
fluid  temperature  rise  is  determined  by  the  Fourier  number  (aT/L2). 

Since  the  maximum  error  for  all  cases  falls  between  Fourier  numbers  of 
0.9  and  2.0,  node  thicknesses  can  be  chosen  sufficiently  small  to  force 
the  Fourier  number  to  a  value  such  that  the  criteria  for  a  long 
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FIG.  17.  Spatial-Truncation  Error  Vs.  Fourier  Number,  Long  Ramp. 
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Fig.  18.  Spatial-Truncation  Error  Vs.  Fourier  Number,  Ultra  Short  Ramp. 
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FIG.  20.  Spatial-Truncation  Error  Vs.  Fourier  Number  for  Ramp 
Function  Temperature  Rise  in  Bounding  Fluid. 
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FIG.  21.  Spatial-Truncation  Error  Vs.  Fourier  Number. 

ramp  exists.  Conversely,  by  choosing  a  sufficiently  large  value  for 
node  thicknesses,  the  Fourier  number  will  be  less  than  0.1  and  the 
conditions  of  the  short-ramp  function  can  be  produced. 

As  previously  noted,  the  analytical  solution  for  the  ramp-function 
problem  shows  that  the  temperatures  after  the  end  of  the  ramp  can  be 
calculated  by  subtracting  a  At  equal  to  the  temperature  rise  from  time 
zero  for  a  time  increment  equal  to  the  elapsed  time  after  the  end  of  the 
ramp  from  the  temperature  generated  if  the  ramp  were  to  continue 
indefinitely.  To  determine  the  temperatures  after  Tj  (the  time 
representing  the  end  of  the  ramp  temperature  rise  of  the  bounding  fluid) , 
successive  temperature  responses  are  generated  as  if  the  ramp  were 
continuing  indefinitely.  However,  to  correct  for  the  fact  that  the 
ramp  did  cease  at  Tj,  the  analytical  solution  subtracts  the  values  of  a 
temperature  response  generated  from  time  zero  by  using  a  time  increment 
equal  to  the  actual  time  increment  after  the  end  of  the  ramp. 

If  this  approach  is  valid  in  the  analytical  solution,  it  should  be 
valid  in  the  equivalent  numerical  solution.  Starting  with  the  numerical 
temperature  calculation  as  a  function  of  time  for  a  given  ramp-function 
problem,  extend  the  generation  of  nodal  temperatures  for  an  indefinitely 
long  ramp  by  merely  „ reversing  the  subtraction  procedure  outlined  in  the 
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analytical  solution.  That  is,  add  the  temperature  rise  found  for  a  time 
increment  after  time  zero  to  the  temperature  for  a  time  equal  to  the 
time  at  the  end  of  the  ramp  plus  the  time  increment.  Thus  an  extension 
of  a  ramp  is  synthesized  by  calculating  the  extensions  from  the  single 
ramp  run.  This  extended  run  was  checked  by  an  independent  numerical 
run  for  the  extended  ramp  and  the  two  results  were  identical.  Thus  with 
a  single  computer  run,  ramps  of  any  length  can  be  synthesized. 

A  further  examination  of  the  analyti cal  solution  reveals  that  the 
final  fluid  temperature  divided  by  the  time  length  of  the  ramp  is  the 
controlling  factor  in  the  calculation  of  nodal  temperatures.  In  other 
words,  the  crucial  variable  in  determining  node  temperatures  is  the 
time  slope  of  the  temperature  rise  in  the  bounding  fluid.  Thus  to 
calculate  temperature  responses  for  other  time  slopes,  the  ratio  of 
the  new  slope  to  the  old  slope  multiplied  by  the  temperature  rise  of 
each  node  in  the  system  would  provide  the  correct  node  temperature  rises 
for  the  new  slope.  For  example,  if  an  analytical  solution  for  all  the 
node  temperatures  caused  by  a  fluid  temperature  rise  of  2500  degrees 
per  second  has  been  calculated,  node  temperatures  for  a  slope  of  250 
degrees  per  second  may  be  determined  by  multiplying  the  node  tempera¬ 
ture  rises  from  the  2500-degree-per-second  results  0.1.  Thus  for  any 
given  problem,  the  node  temperatures  for  the  analytical  solution  may  be 
found  by  multiplying  each  node  temperature  rise  bv  the  ratio  of  the  new 
slope  to  the  old  slope. 

It-  seemed  reasonable  that  if  temperature  responses  in  the  nodes 
could  be  synthesized  for  the  analytical  solutions  by  multiplying 
temperature  responses  by  ratios  of  ramp  slopes,  then  a  similar 
synthesis  method  should  prove  valid  for  the  numerical  solution.  Thus 
a  number  of  runs  for  the  same  geometrical  problem  were  made  such  that 
the  only  variable  was  the  ramp  slope.  The  length  of  the  ramp  in  each 
case  was  held  constant.  As  was  expected  the  numerical  results  from  this 
series  of  runs  were  all  exactly  proportional  to  the  ratio  of  ramp  slopes 
of  temperature  in  the  bounding  fluid. 

By  combining  these  two  separate  synthesis  methods,  it  is  now 
possible  for  a  one-dimensional  equal-node  geometry  and  given  material 
properties  to  synthesize  from  a  single  numerical  run  what  the  node 
numerical  temperature  responses  will  be  for  either  a  longer  or  shorter 
ramp  time  and  for  any  variation  in  ramp  slope. 

It  may  be  noted  that  regardless  of  the  ramp  slope,  the  basic  error 
is  the  same  (Fig.  21).  Whether  or  not  the  ramp  is  25  or  2500  degrees 
per  second,  the  percent  error  is  a  function  of  Fourier  number  only.  It 
should  be  remembered,  however,  that  the  percentage  error  listed  in  Fig. 
21  is  a  percentage  of  the  current  fluid  temperature.  Therefore,  the 
absolute  error  in  the  2500-degree-per-second  ramp  will  be  one  hundred 
times  the  absolute  value  of  the  error  in  the  25-degrees-per-second 
ramp. 
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DELTA-FIRST-DERIVATIVE  CORRELATION 
WITH  SPATIAL-TRUNCATION  ERROR 


The  possibility  of  establishing  a  correlation  between  the  factors 
causing  a  spatial-truncation  error  and  the  error  itself  is  of  interest. 
The  use  of  such  a  correlation  would  allow  the  computer  to  make  a 
numerical  calculation  of  eitner  the  second  derivative  or  a  change  in 
the  first  derivative  to  determine  the  spatial-truncation  error.  This 
approach  might  be  useful  in  practical,  multilayer,  three-dimensional 
problems . 

It  was  first  believed  that  an  evaluation  of  the  second  derivative 
would  be  fruitful.  Upon  making  some  calculations  of  the  second 
derivative  during  the  early  stages  of  a  transient,  it  was  realized  that 
the  second  derivative  was  quite  sensitive  to  the  choice  of  the  thickness 
of  the  individual  nodes.  In  other  words,  as  an  increasingly  thicker 
AX  is  specified,  the  nodes  become  increasingly  less  responsive  to  the 
input  of  energy,  and  the  second  derivative  for  the  same  transient  has 
progressively  smaller  numerical  values  as  AX  is  increased. 

It  was  at  this  point  that  the  decision  was  made  to  see  if  a 
difference  in  the  first  derivatives  might  be  less  sensitive  to  the 
effect  of  increasing  AX.  Calculations  determining  the  magnitude  of  the 
first  derivatives  for  the  first  node  and  the  second  node  seemed  to 
support  the  assumption  that  the  change  in  first  derivatives  might  be 
considerably  less  sensitive  to  node  thickness  than  the  numerical  second 
derivative  under  the  same  conditions.  The  calculation  of  the  first 
derivative  of  the  surface  node  is  made  by  (1)  taking  the  temperature  of 
the  surface,  (2)  subtracting  the  temperature  of  the  first  node,  and 
(3)  dividing  this  difference  by  the  distance  from  the  surface  to  the 
center  of  the  first  node.  The  first  derivative  of  the  second  node  is 
made  by  (1)  taking  the  temperature  of  the  first  node  and  (2)  subtracting 
the  temperature  of  the  second  node,  and  (3)  then  dividing  the  result 
by  the  distance  between  the  first  and  second  node  centers.  The  delta- 
first-derivative  is  the  difference  between  the  two  first  derivatives. 

A  plot  of  this  delta-first-derivative  as  a  function  of  Fourier  number  is 
shown  in  Fig.  22.  The  delta-first-derivative  starts  as  a  large  value 
for  time  zero  and  progressively  reduces  to  smaller  values  as  time 
increases.  Furl 'or  thought  suggested  that  the  spatial-truncation 
error-response  curve  shape,  i.e.,  a  maximum  error  around  a  Fourier 
number  of  1,  could  be  developed  for  the  delta-first-derivative  data. 
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FIG.  22.  Change  in 
First  Derivative  as  a 
Function  of  Fourier 
Number. 
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The  absolute  values  of  the  temperatures  during  the  early  portions  of 
the  transient  are  quite  small.  Thus,  the  effect  of  a  large  change 
in  delta-first-derivatives  would  be  relatively  small  when  multiplied  by 
the  small  temperature  rise.  To  compensate  for  the  large  delta-first- 
derivatives  and  for  the  small  temperature  rises  in  the  early  portion 
of  the  transient,  these  two  factors  were  multiplied  together  to  form  a 
correlation  number.  Figure  23  shows  the  result  of  plotting  this 
correlation  number  versus  the  Fourier  number.  For  comparison  purposes 
the  spatial-truncation  error  expressed  as  a  percentage  of  the  step 
function  in  the  bounding  fluid  is  plotted  as  a  dotted  line.  The 
Fourier  number,  at  which  a  maximum  occurred  in  the  correlation  number, 
coincides  almost  exactly  with  the  maximum  Fourier  number  for  the 
spatial-truncation  error  curve. 

Unfortunately,  at  the  present  time,  there  is  not  a  good  correlation 
in  the  values  of  the  correlation  number  with  respect  to  the  percentage 
values  of  the  spatial-truncation  error  curve.  Presumably,  this  short¬ 
coming  of  the  present  correlation  could  be  overcome  in  a  practical 
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computer  program  wherein  the  correlation  numbers  could  be  generated 
throughout  the  transient.  After  the  computer  had  finished  calculating 
the  entire  transient,  a  subroutine  could  cause  the  computer  to  go  back 
and  match  the  proportional  correlation  number  curve  with  a  spatial- 
truncation  error  curve,  and  to  print  appropriate  temperature  corrections 
at  the  various  time  steps.  This  appears  to  be  a  fruitful  area  for 
further  work  which  could  result  in  a  generalized  error  treatment » 


FIG.  23.  Correlation  Nutter  as  a  Function  cf  Fourier  Number 
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PROGRAMMED  ERROR  PRINTOUT  AND  OPTIMIZATION 


ERROR  PRINTOUT 

The  ONE-D  program  contains  a  section  that  will  compute  the  spatial- 
and  time-truncation  errors.  Sufficient  data  have  not  been  collected  at 
this  time  to  allow  a  similar  error  analysis  for  multilayer  slabs.  Some 
data  are  presented  in  Fig.  11  and  12. 

Since  it  would  take  an  infinite  number  of  error  curves  similar 
to  those  in  Fig.  3  and  4  to  cover  all  possible  nodal  subdivisions 
and  time  increments,  the  error  analysis  can  only  be  approximate. 

Data  from  the  error  curves  for  Biot  numbers  of  “,  10,  3,  1  and 
0.5  were  used  in  a  Chebyshev-polynomial  curve  fitting  program  that 
converted  the  Chebyshev  series  to  its  equivalent  power  series.  The 
power  series  is  of  the  form, 

j  ®  m 

-IV 

j  “  O 

A  10-degree  polynomial  was  used  to  fit  the  spatial  error  curves  and  a 
6-degree  polynomial  for  the  time-truncation  error  curves.  The 
coefficients  for  the  power  series  are  included  with  the  ONE-D  program, 
and  are  listed  in  Appendix  A. 

The  power  series  are  used  to  calculate  the  spatial-truncation  error 
for  a  Fourier  number  between  0.25  and  7.0  except  for  Biot  number  of  0.5, 
where  the  Fourier  number  range  is  0.75  through  7.0.  For  Fourier  numbers 
less  than  0.25  or  0.75,  which  ever  the  case  may  be,  a  value  of  2.5%  for 
the  spatial  error  is  printed.  For  Fourier  numbers  greater  than  7.0,  a 
value  of  -1.3%  for  the  spatial  error  is  printed. 

In  the  case  of  the  time-truncation  error,  the  power  series  are  used 
for  Fourier  numbers  between  0.0  and  10.0.  For  Fourier  numbers  greater 
than  10.0,  the  Fourier  number  is  printed  along  with  a  suggestion  that 
the  time  increment  (DTAU)  be  reduced. 
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The  error  program  tests  the  Biot  number  against  the  Biot  numbers 
for  the  error  curves  used.  Any  Biot  number  falling  between  the 
established  numbers  is  replaced  by  the  next  higher  number,  and  the 
error  is  calculated  on  that  basis. 

The  error  printout  routine  produces  an  error  printout  at  each 
time  step  and  provides  an  option  for  the  user  to  apply  the  error  print¬ 
out  to  ramp  functions.  The  spatial-truncation  errors  in  ramp-function 
temperature  responses  in  bounding  fluids  have  been  Incorporated  in  the 
error  printout.  The  base  error  curve  for  the  long  ramp  has  been 
established  in  terms  of  the  Chebyshev  series  as  a  function  of  Fourier 
number  and  is  stored  in  the  subroutine  for  error  printout.  Thus  in 
rarep-functlon  cases  where  a  long  ramp  is  encountered,  the  user  will 
receive  a  printout  of  the  number  of  degrees  Fahrenheit  due  to  the 
spatial-truncation  error  in  the  first  node.  For  ramps  having  lengths 
shorter  than  the  long  ramp,  that  is,  for  all  ramps  having  Fourier 
numbers  of  less  than  five  at  the  end  of  the  ramp,  the  program  is  so 
devised  that  up  to  the  end  of  the  ramp,  the  program  will  print  errors 
from  the  base  curve.  After  the  end  of  the  ramp,  the  program  will 
automatically  shift  over  to  the  spatial- truncation  error  curve  forming 
the  envelope  for  the  post-ramp  deviations,  as  seen  in  Fig.  21. 

A  FORTRAN  listing  of  the  computer  program  for  error  printout  is 
included  in  Appendix  A. 


AUTOMATIC  TIME-STEP  GENERATOR 

Quite  often  the  user  of  a  transient  heat-transfer  computer  program 
does  not  have  an  accurate  idea  of  a  time  step  to  use  to  minimize  time- 
truncation  errors  at  the  beginning  of  transient.  In  addition,  the  user 
may  not  know  how  many  time  steps  to  specify  to  completely  cover  the 
transient.  In  some  cases,  the  user  can  over  specify  the  number  of  time 
steps  to  be  used  and  waste  computer  time  by  calling  for  printouts  of 
temperatures  after  the  transient  has  been  completed;  or  conversely,  and 
even  worse,  the  user  can  specify  only  a  fraction  of  the  time  steps 
necessary  to  complete  the  transient  and  thus  be  forced  to  go  back  and 
rerun  the  problem,  with  the  resulting  increased  expense  and  delay  before 
obtaining  the  results.  Thus,  if  an  automatic  time-step  generator  were 
available  that  would  keep  the  time-truncation  errors  within  pre¬ 
determined  limits  and  would  cause  the  program  to  cease  generating 
temperature  printouts  when  the  transient  had  completed  a  predetermined 
level  of  response,  a  considerable  savings  in  time  and  improvement  of 
accuracy  could  result.  To  this  end,  an  automatic  time-step  generator 
has  been  developed.  Briefly,  the  basis  for  this  routine  is  that  the 
machine  generates  the  first  time  step  based  upon  the  level  of  accuracy 
the  user  desires,  as  reflected  in  the  Fourier  number.  Once  the  first 
time  step  has  thus  been  determined,  the  program  can  then  proceed  at 
periodic  time  intervals  to  test  each  of  three  nodes  specified:  First, 
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a  node  adjacent  to  the  surface;  second,  a  node  at  an  intermediate  level, 
and  third,  a  node  at  the  deepest  portion  of  the  geometry.  The  test  that 
the  program  runs  is  to  determine  the  temperature  change  for  each  of  the 
three  nodes  during  the  time  step.  If  this  number  of  degrees  is  less 
than  0.5%  of  the  overall  temperature  rise  in  the  bounding  fluid  for  the 
transient,  the  program  will  automatically  double  the  length  of  time  for 
the  succeeding  time  step.  If  the  user  so  desires,  the  percentage  that 
is  used  can  be  changed.  The  final  control  of  the  program  will  permit 
the  user  to  specify  the  terminal  temperature  in  any  node  desired. 


NWCTP  5143 


ONE-D  ACCELERATION  OF  THREE-DIMENSIONAL 
TRANSIENT  PROBLEMS 


The  more  complex  three-dimensional  heat-transfer  problems  involv¬ 
ing  very  slow  convergence  pose  the  problems  of  excessive  machine  time 
and  possible  inaccuracies  due  to  the  convergence  characteristics  of 
the  iterative  processes  used.  These  iterative  processes  depend  upon 
some  prescribed  limit  for  termination,  and  the  value  of  this  iteration 
limit  for  maximizing  accuracy  and  for  minimizing  machine  time  is  not 
known  a  priori.  To  reduce  machine  time  and  to  have  less  dependence 
upon  the  sensitivity  of  the  iteration  limit,  a  chain  of  nodes  could  be 
established  in  a  psuedo  ONE-D  pattern  extending  from  the  surface  of  the 
three-dimensional  geometry  to  the  deepest  portion  of  the  geometry. 

These  chains  of  nodes  could  be  solved  directly  and  could  establish  a 
temperature  field  as  the  beginning  point  for  a  computer  program  for 
three-dimensional  thermal  response,  THT-B.  This  iterative  process 
would  then  be  reduced  to  essentially  lateral  heat  transfer  and,  in  many 
cases,  be  capable  of  convergence  to  the  final  solution  within  a  few 
iterations.  Preliminary  work  on  this  approach  consisted  of  establishing 
a  very  slowly  converging  two-dimensional  model.  Temperatures  for  all 
nodes  are  evaluated  at  the  beginning  of  each  time  step  by  the  separate 
application  of  the  one-dimensional  TRIDAG  technique  to  two  columns  of 
nodes.  The  entire  temperature  field  of  the  model  is  then  solved  by 
using  the  Gauss-Seidel  method.  The  results  are  summarized  in  Table  2. 
These  preliminary  results  indicate  that  it  may  be  profitable,  from  the 
dual  standpoint  of  accuracy  and  reduction  in  machine  time,  to  establish 
a  subroutine  wherein  the  user  may  elect  in  specific  three-dimensional 
problems  to  establish  one-dimensional  chains  of  nodes  for  which  a 
tridiagonal  solution  would  automatically  be  used  at  the  beginning  of 
each  time  step  to  establish  the  approximate  temperature  field  before  an 
iterative  technique  is  utilized. 
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TABLE  2.  Comparison  of  Results  of  the  THT-B  Computer  Program 
With  Those  of  the  THT-B  and  the  ONE-D  Programs  Combined. 
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SUMMARY  AND  CONCLUSIONS 


To  summarize  the  work  described  in  this  report,  all  of  the 
component  projects  were  directed  (1)  to  improve  the  accuracy  of  the 
transient  numberical  solutions  currently  available,  and  (2)  to  reduce 
the  machine  time  necessary  for  obtaining  accurate  transient-temperature 
histories. 

This  study,  consisting  of  various  projects  devcted  to  different 
aspects  of  numerical  solutions  of  transient  heat-transfer  problems,  has 
arrived  at  an  evaluation  of  spatial-  and  time- truncation  errors  for 
homogeneous  slabs  using  uniform  node  spacing  and  being  exposeu  to 
bounding  fluids  that  have  either  a  step-function  or  a  ramp-function 
change  in  temperature.  The  errors  involved  in  thin- thick  geometries 
of  a  thin  layer  of  highly  conductive  material  in  contact  with  a  thick 
layer  of  insulating  material  are  presented.  With  the  information 
provided  by  the  curves,  a  user  can  ensure  that  the  accuracy  of  his 
transient  temperature  analysis  and  design  will  fall  within  any  pre¬ 
determined  level  of  accuracy.  Whenever  possible,  the  curves  reflect 
the  worst  possible  case;  therefore,  in  practical  problems  where 
conditions  not  as  severe  as  the  study  conditions  occur,  a  conservative 
evaluation  of  the  errors  will  automatically  result.  With  respect  to  the 
speeding-up  of  the  solutions,  particularly  with  regard  to  computer  time 
used,  several  acceleration  techniques  are  presented  with  some  evaluation 
as  to  their  effectiveness.  In  addition,  a  direct  solution  applicable 
to  one-dimensional  problems  is  provided  that  reduces  computer  solution 
time  when  compared  to  the  usual  iterative  solutions.  As  a  fringe 
benefit  of  this  study,  a  number  of  the  small  computer  programs  used  in 
the  error  study,  acceleration  and  other  portions  of  the  study,  have  been 
combined  into  one  generalized  transient  heat-transfer  one-dimensional 
program  titled,  ONE-D.  This  generalized  program  should  prove  to  be  of 
considerable  help  to  the  heat-transfer  transient  designer,  for  a 
number  of  options  are  available. 

By  examining  the  error  response  of  ramp  functions,  the  user  now  has 
a  mathematical  model  that  is  considerably  closer  to  practical  problems 
than  is  the  step-function  model  previously  used.  The  error  printout 
response  both  to  the  step  function  and  to  the  ramp  function  at  the 
choice  of  the  user.  The  principle  of  utilizing  the  ONE-D  solution  for 
improving  both  the  accuracy  and  the  machine  time  aspects  of  slowly 
converging  three-dimensional  programs  will  be  beneficial  to  the  user. 
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In  like  manner,  the  automatic  time-step  generator  will  also  afford  the 
user  convenience,  accuracy  and  reduced  machine  time.  The  correlation 
of  the  delta-first-derivative  with  spatial-truncation  error  may  prove 
to  be  the  basis  for  a  generalized  error  treatment  in  which  the  mathe¬ 
matical  model  will  be  the  problem  reduced  to  numerical  terms  instead 
of  modeling  real  problems  with  the  one-dimensional  homogeneous  equal 
AX  geometry. 

In  conclusion,  the  user  is  now  in  a  position  (1)  to  define  much 
more  accurately  the  areas  of  real  transient  heat-transfer  problems  in 
which  error  may  be  considerable;  (2)  to  change  his  input  data  to  make 
sure  that* his  results  are  reasonably  accurate;  and  (3)  to  make 
maximum  effective  use  of  the  transient  heat-transfer  programs  available 
utilizing  the  techniques  outlined  here  such  as  the  ONE-D  acceleration 
and  the  automatic  time-step  generator. 
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Section  1 

GENERAL  HEAT-TRANSFER  PROGRAM  “ONE-D”  AS  APPLIED 
TO  STEP-FUNCTION  BOUNDARY  CONDITIONS 


A  number  of  specialized  small  computer  programs  were  developed  in 
the  course  of  this  study  to  investigate  various  aspects  of  the  step- 
function  boundary  conditions.  These  small  programs,  e,g.,  error 
printout,  iteration-acceleration  method,  direct  solutions,  spatial- 
and  time-truncation  error,  and  thin-thick  geometries,  were  combined 
into  a  single  computer  program  called  "ONE-D",  coded  in  FORTRAN  for 
the  IBM  1620  computer. 

This  program  is  capable  of  handling  a  one-dimensional  problem  of 
fifty  nodes;  however,  the  number  of  nodes  may  be  increased  merely  by 
enanging  the  DIMENSION  statement.  The  equations  in  the  program  are  set 
up  such  that  one  surface  of  the  slab  is  exposed  to  a  fluid  while  the 
ocher  surface  of  the  slab  is  insulated. 

The  program  contains  a  number  of  options  that  may  be  called  by 
the  user.  As  an  example,  the  program  can  be  used  for  a  composite  slab, 
made  up  of  three  materials,  as  well  as  for  a  homogeneous  slab.  The 
boundary  conditions  of  fluid  temperature  distribution  of  the  slab  can 
be  included  in  the  input  or  can  be  considered  constant.  The  solution 
of  the  eciuations  in  the  program  can  be  accomplished  in  a  number  of  ways 
such  as  Gauss-Seidel,  accelerated  Gauss-Seidel,  or  TRIDAG,  a  direct 
solution  for  tridiagonal  matrices. 


CONTROL  VARIABLES 

The  following  variable  names  are  used  for  the  various  options: 

KODE  The  variable  name  KODE  is  used  to  determine  if  the  initial 
temperature,  of  the  body  is  constant  or  variable.  If  KODE 
is  1,  the  initial  temperature  is  constant.  If  KODE  is  not  1, 
then  the  initial  temperature  must  be  specified  for  each  node. 
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K0DE1 


The  variable  name  K0DE1  is  used  to  determine  if  the  fluid 
temperature  and  film  coefficient  are  constant  or  variable. 

If  K0DE1  is  1,  the  fluid  temperature  and  film  coefficient 
are  constant.  If  K0DE1  is  not  1,  then  the  fluid  temperature 
and  film  coefficient  for  21  values  of  time  must  be  specified 
in  the  input. 

The  variable  name  SEL1  is  used  to  determine  if  equations  are 
solved  by  an  iterative  method  or  by  Tridag.  If  SRL1  is  1, 
the  solution  is  iterative.  If  SEL1  is  not  1,  then  Tridag 
is  used. 

The  variable  name  SEL2  is  used  to  determine  the  type  if 
iterative  solution  to  be  used.  If  SEI.2  is  1,  then  regular 
Gauss-Seidel  iteration  is  used.  If  SEL2  is  not  1,  then 
accelerated  Gauss-Seidel  is  used. 


A  comparison  of  N0DES1  and  NODES  determines  whether  the  slab  is 
homogeneous  or  multilayer.  If  N0DES1  equals  NODES,  the  slab  is  homo¬ 
geneous,  and  the  properties  only  of  the  first  layer  are  read.  Node¬ 
spacing  within  a  layer  must  be  equal,  but  the  node-spacing  of  the 
different  layers  need  not  be  the  same.  This  allows  for  small  node¬ 
spacing  in  thin  layers  and  larger  node-spacing  in  thick  layers. 

Provisions  have  been  made  to  allow  for  contact  resistance  between 
layers  or  between  nodes.  The  variable  ZH  is  used  for  this  purpose. 

ZH1  is  used  as  the  conductance  between  the  nodes  of  the  three  layers. 
ZH2  is  the  conductance  between  layer  1  and  layer  2,  while  ZH3  is  the 
conductance  between  layer  2  and  layer  3. 

The  other  variable  names  are  defined  at  the  beginning  of  the 
program  listing. 

Two  matrix  methods,  iterative  and  tridiagonal  are  used  in  this 
study  to  solve  the  equations  formed  by  the  heat  balance  for  each  node 
in  the  slab. 

The  program  contains  two  sections  that  use  the  iterative  method. 
Section  1  uses  regular  Gauss-Seidel  iteration  to  solve  the  set  of 
simultaneous  equations.  Section  2  uses  Gauss-Seidel  iteration  with 
either  or  both  acceleration  techniques,  Wegstein  and  constant 
acceleration  factor.  The  respective  sections  in  the  program  listing 
are  noted  with  comment  cards. 

The  Wegstein  technique  causes  the  programming  to  be  more  complex 
than  either  TRIDAG  or  regular  Gauss-Seidel  since  the  temperatures  for 
all  nodes  for  the  latest  three  iterations  must  be  stored.  In  the 
Wegstein  technique  the  following  general  equations  are  used: 
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where 


T11  =  Tn  +  Q  (Tn_i  -  Tn) 


T  =  the  accelerated  temperature  for  the  nth  iteration  of  any 
node 

Tn  =  the  temperature  for  the  nth  iteration 
n_  ^ 

T  =  the  temperature  for  the  previous  iteration,  and 

t"  -  t"'1 

A  =  - : - ^ 

Tn-1  _  Tn-2 


Values  of  Q  are  negative  since  values  of  A  are  restricted  to  0<A<1  to 
insure  only  positive  values  of  acceleration.  As  A  approaches  1,  the 
value  of  Q  can  become  quite  large.  Therefore  Q  is  restricted  to  a 
maximum  value  of  100. 

The  Wegstein  technique  is  applied  to  each  node,  so  individual 
values  of  A  and  Q  must  be  calculated.  The  alternative  is  to  apply  a 
constant  acceleration  factor  to  all  nodes;  this  is  done  by  using  the 
variable  FACTOR.  The  variable  FACTOR  is  used  in  the  same  manner  as  is 
Q;  therefore,  it  should  be  negative.  Note  that  if  FACTOR  is  zero, 
there  is  no  acceleration;  thus  regular  Gauss-Seidel  results.  The 
regular  Gauss-Seidel  section  of  the  ONE-D  program  could  have  been 
eliminated  and  the  FACTOR  routine  substituted  by  setting  FACTOR  =  0, 
but  it  was  more  convenient  to  include  it  in  the  program  when  accelera¬ 
tion  study  runs  were  made. 

In  the  program  a  constant  acceleration  factor,  FACTOR,  is  applied 
to  all  nodes  until  ISTOF  is  reached.  NOGS  iterations  of  Gauss-Seidel 
are  made  until  IAPPLY  is  reached;  at  which  iteration,  the  Wegstein 
technique  is  applied.  The  Wegstein  technique  is  repeated  after  the 
interval  INTER. 


If  the  Wegstein  technique  is  applied  too  soon  or  too  often,  the 
solution  may  be  slowed  down.  Past  experience  has  shown  that  the 
initial  application  should  be  made  on  about  the  iteration  which  is 
proportionate  in  number  to  the  depth  of  the  deepest  node  from  the 
boundary.  For  example,  in  the  20-node  problem,  the  initial  application 
of  the  Wegstein  technique  should  be  on  about  the  tenth  iteration.  The 
criterion  suggested  for  frequency  of  application  of  the  technique  is 
given  as  about  one-half  the  number  of  iterations  used  for  the  initial 
application  of  acceleration. 
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By  using  proper  values  of  the  variables  FACTOR,  ISTOP,  NOGS,  and 
INTER,  it  is  possible  to  reduce  the  number  of  iterations  compared  to 
Gauss-Seidel  by  as  much  as  63%.  Obviously  there  is  a  large  number  of 
combinations  possible,  some  of  which  are  undesirable.  If  a  constant 
acceleration  factor  is.  to  be  used  without  Wegstein  acceleration,  the 
literature  indicates  that  FACTOR  should  be  between  0  and  -1,  possibly 
-0.5  or  -0.6.  When  used  with  Wegstein,  FACTOR  may  be  a  larger  numerical 
number  but  still  negative.  From  a  series  of  twelve  computer  runs  made 
on  the  sample  20-node  problem,  it  was  found  Chat  values  from  -11.4  to 
-1.5  for  FACTOR  gave  the  minimum  number  of  iterations. 

Results  of  the  20-node  problem  with  constant  initial  temperature, 
constant  fluid. temperature,  and  film  coefficient  subjected  to  a  step* 
change  in  the  boundary  fluid  are  as  follows : 


Gauss-Seidel  Iteration  Technique 


FACTOR3 

h 

J 

Number  of 

ISTOP 

N0GSC 

INTER 

iterations 

0 

60 

lit 

•  •  • 

46 

-1 

60 

•  •  * 

«  •  % 

35 

-1 

3 

0 

3 

42 

-1 

5 

0 

5 

24 

-1 

6 

0 

3 

40 

-1 

10 

3 

6 

24 

-1.3 

10 

3 

6 

17 

-1.4 

10 

3 

6 

17 

-1.5 

10 

3 

6 

17 

-1.5 

10 

0 

6 

20 

-1.6 

10 

3 

6 

33 

0 

10 

0 

6 

27 

-0.5 

50 

. . . 

31 

3 

Acceleration  factor 

for  all  nodes 

in  system. 

Number 

of 

.iterations 

that  FACTOR 

is  to  be  applied. 

CNumber 

of 

iterations 

between  application  of 

unaccelerated 

Gauss-Seidel 

iterations 

and  Wegstein 

accelerations . 

Number 

of 

Gauss-Seidel  iterations 

between  applications  of 

Wegstein  accelerations. 

The  TRIDAG  solution  needs 

no  special  instructions 

or  data  cards 

The  solution  is  rapid  and  contains  no  convergence  errors  but  is  limited 
to  linear  equations. 
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ARRANGEMENT  OF  INPUT 


Data  Card  1: 

FORMAT  (514,  5F10.3),  NODESl,  NODES 2,  NODES,  K0DE,  KODE1 , 
SEL1,  SEL2,  DTAU,  ZNUM,  ZLIM.  This  card  must  be  used  for 
all  options. 


Data  Card  2: 


FORMAT  (F10.0,  315),  FACTOR,  ISTOP,  NOGS,  INTER.  This  card 
is  read  only  if  accelerated  Gauss-Seidel  is  to  be  used  (SEL2  =  1) . 

Data  Card  3: 


FORMAT  (4F10.3,  E10.3),  DEI.X1,  ZK1,  ZRH01,  ZC1,  ZH1.  This 
card  must  be  used  for  all  options. 

Data  Card  4: 


FORMAT  (4F10.3,  E10.3),  DELX2 ,  ZK2,  ZRH02,  ZC2,  ZH2 .  This 
card  is  read  if  slab  is  made  up-  of  two  or  three  lavers  (N0DES7, 
NODESl) . 

Data  Card  5: 


FORMAT  (4F10.3,  E10.3),  DELX3 ,  ZK3 ,  ZRH03,  ZC3 , 
card  is  read  if  slab  is  made  up  of  three  layers 


ZH3 .  This 
(NODES 7 ,  NODES 2) . 


Data  Card  6-8: 


FORMAT  (8F10.3),  TS .  These  cards  are  used  for  reading  variable 
fluid  temperature  (K0DE1  =  1) .  Twenty-one  values  of  TS  are 
read.  These  values  are  temperatures  of  the  fluid  for  20  equal 
time  steps  covering  the  time  range  for  the  transient. 


Data  Card  9-11: 


FORMAT  (8F10.3),  HFVAR.  These  cards  are  used  for  reading 
variable  film  coefficient  (K0DE1  =  1) .  Twenty-one  values  of 
HFVAR  are  read.  These  values  are  film  coefficients  for  20  equal 
time  steps  covering  the  time  range  for  the  transient. 

Data  Card  12: 

FORMAT  (2F10.3)  TC,  HFC.  This  card  is  read  if  the  fluid 
temperature  and  the  film  coefficient  are  constant  (K0DE1  ~  1). 
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Data  Card  13-19 : 

FORMAT  (8F3.0.3),  TPR.  These  cards  are  used  to  read  the 
nodal  temperatures  if  the  slab  has  an  initial  temperature 
distribution  (KODE  =  1).  The  number  of  data  cards  is 
dependent  on  the  number  of  nodes  in  the  system  with  a 
maximum  of  50  nodes ; 

Data  Card  20: 

FORMAT  (F10.3)  TAU.  This  card  is  read  to  indicate  time  for 
variable  initial  temperature  distribution  (KODE  =  1). 

Data  Card  21: 

FORMAT  (F10.3)  TIN.  This  card  is  read  when  slab  has  a 
constant  initial  temperature  (KODE  =  1). 

Data  Card  22: 

FORMAT  (15),  ITERS,  This  card  is  read  when  the  problem  is 
solved  by  an  iterative  method  (SEL1  =  1) .  This  puts  an 
upper  limi on  the  number  of  iterations  allowed  to  converge. 

Data  Card  23-47 : 

These  cards  are  supplied  with  the  program.  Their  purpose  is 
to  supply  the  coefficients  of  the  power  series  that  describe 
the  spatial  and  time  truncation  error  curves .  These  cards  are 
read  only  if  the  slab  is  homogeneous  (NODES 1  =  NODES)  and  for 
the  first  time  step  (TAU  =  DTAU) . 
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C  ONE -DIMENSIONAL  HEAT  TRANSFER  PROGRAM,  ONE-D,  CODED  BY  VAN  TASSEL 

C  AND  STANLEY 

C 

C  AR.BR.CR.DR  *  COEFFICIENT  ARRAYS  CONTAINING  THE  SUB' DIAGONAL*  DIA- 
C  GONAL,  SUPER-DIAGONAL,  AND  RIGHT  HAND  ELEMENTS  OF 

C  THE  TRIDIAGONAL  SYSTEM 

C  BIOT  *  BtOT  NUMBER 

C  COEFF  *  COEFFICIENTS  OF  POWER  SERIES  USED  IN  ERROR  ANALYSIS 
C  DELX  »  SPATIAL  INCREMENT,  INCHES 

C  DTAU  *  TIME  STEP.  SECONDS 

C  FACTOR  *  CONSTANT  OVER  RELAXATION  FACTOR (ZERO  FOR  REGULAR  GAUSS' 

C  SEIDEL) . 

C  FO  «  FOURIER  NUMBER 

C  HFC  »  CONSTANT  FILM  COEFFICIENT,  BTU/HR-SO  FT*F 
C  HFILM  *  FILM  COEFFICIENT  AT  SPECIFIC  TIME  UNDER  CONSIDERATION, 

C  BTU/HR-SO  FT-F 

C  HFVAR  *  TIME  VARIA8LE  FILM  COEFFICIENT,  BTU/HR-SO  FT-F 

C  INTER  •  NUMBER  OF  REGULAR  GAUSS-SESDEL  ITERATIONS  BETWEEN  SUCCES- 

C  SiVE  APPLICATIONS  OF  THE  WEGSTEIN  ACCELERATION. 

C  I  STOP  *  NUMBER  OF  ITERATIONS  THAT  THE  CONSTANT  OVER  RELAXATION 
C  FACTOR  IS  APPLIED. 

C  ITERS  «  ALLOWABLE  NUMBER  OF  ITERATIONS 

C  KSWPS  *  NUMBER  OF  ITERATIONS  PER  TIME  STEP 

C  NODES  I  *  NUMBER  OF  NODES  IN  FIRST  LAYER 

C  N0DES2  *  NUMBER  OF  NODES  IN  SECOND  LAYER 

C  NODES  ■  TOTAL  NUMBER  OF  NODES  IN  THE  SYSTEM 

C  NOGS  «  NUMBER  OF  REGULAR  GAUSS-SEIDEL  ITERATIONS  APPLIED  BEFORE 
C  WEGSTEIN  ACCELERATION  IS  APPLIED. 

C  TAU  «  TIME,  SECONDS 

C  T  •  TEMPERATURE  OF  NODE  AT  BEGINNING  OF  ITERATION,  F 

C  TC  *  CONSTANT  FLUID  TEMPERATURE,  F 

C  TF  *  FLUID  TEMPERATURE  AT  SPECIFIC  TIME  UNDER  CONSIDERATION,  F 

C  TIN  •  CONSTANT  INITIAL  SLAB  TEMPERATURE,  F 

C  TNEW  *  TEMPERATURE  OF  NODE  AT  END  OF  ITERATION  OR  TIME  STEP,  F 
C  TPR  «  TEMPERATURE  OF  NODE  AT  PREVIOUS  TIME  STEP,  F 
C  TS  *  TIME  VARIABLE  FLUID  TEMPERATURE,  F 

c  7 r  .  cpcrtFtr  up.r  DTII/I  RM-P 

C  ZH  *  CONDUCTANCE  BETWEEN  NODES, TO  ALLOW  FOR  CONTACT  RESISTANCE 
C  ZK  «  THERMAL  CONDUCTIVITY,  BTU/HR-FT-f 

C  ZRHO  «  DENSITY,  LBM/CU.FT, 

C  ZLIM  *  ITERATION  LIMIT 

C  ZNUM  «  NUMBER  OF  TIME  STEPS  DESIRED 

C 

DIMENSION  TPR 1 50 1 , T ( 50 1 , TNEW (50) ,TT( 50 ,3 ) ,FTl 50 ,3 ) , TOLD (SO, 3 1 , 

I  C(50 ) ,Dl 50 ) , AR( 50 ) ,BR(50 ) ,CR(50 I ,DR(S0 ) , BETA (50 ) .GAMMA (50 ) ,ZK<  50 
2) ,ZC(50 ) , ZRHO (50 ! , DELX (50 ) ,ZH(50 ) ,TS(2t ) ,HFVAR(2I ) ,COEFF(5, 1 1 » 

10  FORMAT ( 5 1 4 , 5F 1 0 . 2 ) 

20  FORMAT ( IH  ,UX,7HFILM  H*,FI0.2,7X, I IHBOUND  TEMP* ,F8.2,2H  F) 

30  FORMAT  (8F10.3) 

40  FOriAT(2FI0.3) 

50  FORMAT ( IH0.2XMHNODE,  I IX.7HDELTA  X , MX , I HK , I 6X , 7HDENS I TY , I 5X , IHC, 
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! 

i 


i 

i 


i 


1 ISX.2HHC) 

80  FORMAT ( IH  ,  IX,  13,5X,FI4.3,SX,F!4,3,5X,FI4.3,5X,FI4.3,9X,EI0.3) 

70  FORMAT!  IH  .  I 1XY4HTIME.FI2.3, 10H  SECONDS  , / ,  I IX, I  IHITERATIONS*  ,13) 
80  FORMAT  I  IH  .SFH .3,SX,SFI 1.3) 

90  FORMAT! IH  , 18X, I3,20X,FI0,5) 

100  FORMAT! i HO, 24H ITERATION  LIMIT  EXCEEDED) 

110  FORMAT! SH0, I IX, 1 4HINI TIAL  TIME*  ,FI0.0,i0H  SECONDS) 

1.20  FORMAT!  IH0.78HFILM  COEFFICIENT, BEGINNING  AT  TIME  ZERO, WITH  TIME  IN 
ICREM£NT*TOTAL  TIME  SPAN/20) 

130  FORMAT! 1 HO. 73HFLUID  TEMP. .BEGINNING  AT  TIME  ZERO, WITH  TIME  INCREME 
1  NT* TOTAL  TIME  SPAN/20) 

140  FORMAT ( F i 0 . 0 , 3 1 5 ) 

ISO  FORMAT! IS) 

180  FORMAT (8E 1 0.3) 

170  FORMAT «5E IS. 8) 

180  FORMAT ( 4F 1 0  <  3 , E 1 0 . 3 ) 

190  FORMAT! IH0.49HMAXIMUM  SPATIAL  TRUNCATION  ERROR  IS  *2.5  PER  CENT) 
200  FORMAT! IH0.37HMAXIMUM  SPATIAL  TRUNCATION  ERROR  IS  -.F6.I.9H  PER  CE 
Tl 

210  FORMAT! I HO, 77HT I ME  STEP  IS  TOO  BIG,  SUGGEST  USING  SMALLER  FOURIER 
UMBER.  FOURIER  NUMBER  IS, FI 0.2) 

220  FORMAT! IH0,34HMAX1MUM  TIME  TRUNCATION  ERROR  IS  -.F6.I.9H  PER  CENT) 

230  FORMAT! IH0.49HMAXIMUM  SPATIAL  TRUNCATION  ERROR  IS  -1.3  PER  CENT) 
240  FORMAT! I H  , 1IX.IIHNODE  NUMBER, I8X, I 7HTEMPERATURE  DEG  F) 

250  FORMAT! I  HO, I9HFILH  COEFF. ‘CONST. « ,FI 0.3 ) 

260  FORMAT! !H0, ISHFLUID  TEMP. *CONST. * ,FI 0.3 ) 

270  FORMAT! IHO, I IX.4HTIME.FI2.3, !0H  SECONDS) 

280  FORMAT! FI 0,3) 

READ  1 0 . NODES  I . NODES2 , NODES , MODE , KOBE  I , SELI , SEL2, DTAU . ZNUM . ZL I M 
IF! SEL2- I • J242.243.242 

C  IF  SEL2  IS  I,  SOLUTION  BY  ORDINARY  GAUSS-SEIDEL,  IF  NOT,  SOLUTION 
C  BY  ACCELERATED  GAUSS -SEIDEL, FAC TOR, I  STOP, NOGS. AND  INTER  MUST  BE 
C  READ. 

242  READ  1 40, FACTOR. I STOP, NOGS. INTER 

243  PRINT  50 

C  READ  PROPERTIES  OF  MATERIALS,  MAXIMUM  OF  THREE  MATERIALS 
READ  I80.DELXI ,ZKI .ZRHOI ,ZC1 ,ZHI 
IF  INODES  I -NODES  123 1. 232, 46 

C  DETERMINE  IF  HOMOGENEOUS  MATERIAL,  IF  NODESI ‘NODES  -«  HOMOGENEOUS. 

231  READ  180, DELX2 , ZK2 , ZRH02 , ZC2 , ZH2 
I F ( NODES2 -NODES ) 233 ,232,46 

233  READ  l80.DELX3,ZK3,ZRHO3,ZC3,ZH3 

C  SET  PROPERTIES  OF  EACH  NODE  IN  EACH  LAYER 

232  DO  93  N* I, NODES 
DELX(N) *DELXI 
ZK!N)«ZKI 
ZRHOI N) * ZRHOI 
ZC(N) *ZCI 

93  ZHIN) *ZHI 

I F ( NODES  I -NODES ) 234 . 236 , 46 

234  ZH I  NODES  I )*ZH2 
1 ‘NODES  I ♦ I 

DO  102  N* I .NODES2 
DELXIN) ‘DELX2 
ZKIN) *ZK2 
ZRHOI N) ‘ZRH02 


4 
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102 

235 


103 

236 

4 


C 

C 

C 

C 

C 


237 


238 


239 

241 

S 


C 

C 

I 


2 


3 


7 

33 

C 

C 

C 

36 

6 


ZC<N)*ZC2 
ZHlN) *ZHI 

1 F i N0DHS2 -NODES  1235,236,46 

ZHtNODES21*ZH3 

I *NODES2* 1 

DO  103  N* l, NODES 

DELX(N) *DELX3 

ZK(N)*ZK3 

ZRHOlM) *ZRH03 

ZC(N) *ZC3 

ZHlN) *ZHI 

DO  4  N» I, NODES 

C (N) *300, *ZC(N) ‘ZRHOIN1 *DELX(N) /DTAU 
PRINT  60,N,DELX(N) , ZKl N) , ZRHOi N) ,2C(N) ,ZH(N) 

IFlKODEI-l 1237,238,237 

IF  NODE  I  IS  I,  THE  FLUID  TEMPERATURE  AND  FILM  COEFFICIENT  ARE  CON¬ 
STANT.  IF  NODE  IS  NOT  I,  VARIABLE  FLUID  TEMPERATURES  AND  FILM 
COEFFICIENTS  ARE  READ; 

READ  AND  PRINT  VARIABLE  FLUID  TEMPERATURE  AND  FILM  COEFFICIENT 
READ  30, (TS(I), 1*1,211 
PRINT  130 

PRINT  80, (TS( 1 1,1*1,211 
READ  30, IHFVARI 11,1*1,21) 

PRINT  120 

PRINT  8C , IHFVARI 11,1*1,2!) 

GO  TO  241 
READ  40.TC.HFC 
DO  239  1*1,21 
TS( 1 1 *TC 
HFVARl 1 1 »HFC 
NP* NODES -I 
DO  5  N.*I,NP 

DlN) *  I ./< DELXlN) /! 24. *ZK<N) 1 *DELXIN* 1 1/(24, *ZK(N* I  1 1  *  I ,/ZHlN) ) 
IFIKODE-1 1 1 ,3, 1 

IF  KODE  IS  ! ,  ALL  NODES  ARE  SAME  TEMPERATURE.  IF  MODE  IS  NOT  EQUAL 
TO  I  NODAL  TEMPERATURES  ARE  READ  IN. 

DO  2  N*I, NODES 
READ  30.TPRIN) 

PRINT  80.TPRIN) 

TIN) *TPR(N1 
READ  280, TAU 
PRINT  110, TAU 
60  T(  33 
READ  280. TIN 
DO  7  N* I o NODES 
TPRlN)  *  TIN 
TIN) *TPRIN) 

TAU*0.0 

IFISELI -  1 . 135,36,35 

IF  SELI  IS  I,  THE  SOLUTION  IS  BY  EITHER  GAUSS-SEIDEL  OR  BY  ACCEL¬ 
ERATED  GAUSS -SEIDEL  AND  THE  MAXIMUM  ALLOWABLE  NUMBER  OF  ITERATIONS 
MUST  BE  READ.  IF  SELI  IS  NOT  I,  SOLUTION  IS  BY  TRIDAG. 

READ  150, ITERS 

TAU*TAU»DTAU 

KSWPS«0 


o: 


l 

I 
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o  o 


MOUNT  *  0 
U«l 

SK* I • *TAU*20»/( ZNUK*DTAU » 

M*SM 

RM*M 

I F I TAU - ZNUM  *DTAU >  254 , 255 , 46 

C  DETERMINE  FILM  COEFFICIENT  AND  FLUID  TEMPERATURE  FOR  TIME  UNDER 
C  CONSIDERATION. 

284  TF*TS(K) ♦ ITSIM* I )-?SlK) ) *(SK-RK) 

HF I LM*HF VAR ( K I ♦ ( HFVAR I  M* I ) -HFVAR < K  )l  *  I SK-RK I 
GO  TO  256 
2SS  TF*TSIM) 

HF1LM*HFVAR(K) 

2E8  A*  I ./-( It . /HFILMJ  *DELX(  1 1/124.  *ZK(  1 1  >8 
PRINT  20.HFILM.TF 
IFISEL2-I • 16) ,8,61 
C 

C  REGULAR  GAUSS-SEIDEL  ITERATIVE  SOLUTION 

8  IFIKSKPS- ITERS )9»9, 28 
3  DO  28  N* I .NODES 
IF (N* 1 1 12, 12, I  I 

11  IF(N-NODES) 14,13,13 

12  TNEKl I ) *  I A«TF*CC I J*TPR( I )*Dt I )*T<2> )/(A*C( 1 1 »0( t  > I 
GO  TO  IS 

13  TNEHI NODES) *( C (NODES) • TPRI NODES ) *D I  NODES- 1 1 *T I  NODES- I  I J/lClNODES) ♦ 
ID( NODES- I ) ) 

GO  TO  IS 

14  TNEKl N) *  1C (N) *TPR(N) *D(N- I ) ‘TIN- 1 ) *DlN) *TCN* I ) )/(C(N) *D(N- 1 ) *D(N) ) 

15  TEMP*TNEK(N) -TIN) 

IF(TEMP) 16,17,1? 

16  TEMP* (-TEMP) 

I?  IF(TEMP-ZLIM) 18,13,19 

18  KOUNT* MOUNT* I 

I F ( MOUNT -NODES ) 25 , 26 , 28 

19  MOUNT *0 

25  TIN) *TNEH(N) 

MSWPS*MSWPS* I 
GO  TO  8 

C  END  OF  GAUSS-SE10EL  SOLUTION 

26  TIN) *TNEK(N1 
24  PRINT  70,  TAU.MSHPS 

PRINT  240 

PRINT  90,  IN,  TIN), N* I , NODES) 

21  DO  27  N* I, NODES 

27  TPRCN) -TIN) 

IFI NODES  I -NODES) 47, 198,46 

28  PRINT  100 
CALL  EXIT 

SPATIAL  AND  TIME  TRUNCATION  ERROR  ANALYSIS 

198  IF(TAU*DTAU)227,227,226 

C  READ  COEFFICIENTS  OF  POKER  SERIES  FOR  SPATIAL  TRUNCATION  ERROR 
227  DO  199  M* I ,5 

199  READ  170,  (COEFFIM, I ) , l* I , I  I  ) 

BIOT  *  HFlLM*DELXI/24./ZMI 

FO  »  ZMI *TAU/(6.25*ZCI *ZRHOI ‘DELXI **2) 
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c 


c 


c 


c 

c 
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DETERM 1NE  SPECIAL  CONDITIONS 
IFIFO-7. 1202,202,201 

201  PRINT  230 
GO  TO  225 

202  1FIBIOT-. 5)203, 203, 205 

203  IF1FO-. 75120A, 215,216 

204  PRINT  190 
GO  TO  225 

205  IFIFO-. 25)206, 207, 207 

206  PRINT  190 
GO  TO  22S 

DETERMINE  WHICH  CURVE  TO  USE,  BASED  ON  BIOT  NUMBER  ' 

207  IFIBIOT-IO. )2! I ,209,206 

208  K  *  5 

GO  TO  217 

209  K  *  4 

GO  TO  217 

21!  IFiBIOT-3. 1213,212,209 

212  K  *  3 

GO  TO  217 

213  1FIBIOT *1.1215,214,212 

214  K  *  2 

GO  TO  217 

215  IFIBIOT-. 5)216, 216, 214 

216  K  *  I 

217  SUM  *  COEFFIK, I ) 

IF 1 1  *  I  1 ) 223 , 2 18,218 

218  DO  219  1*1,10 

219  SUM  «  SUM » COEFFIK, I ♦ U  »FO»*I 
PRINT  200  .SUM 

READ  COEFFICIENTS  OF  POWER  SERIES  FOR  TIME  TRUNCATION  ERROR 

225  DO  221  K*  I  ,5 

221  READ  170,  I COEFFIK, 11,1*1,7) 

FO  *  ZKI *DTAU/(6.25*ZCI *ZRHOI *DELXI »»2) 

IFIFO-IO. 1207,207,222 

222  PRINT  210,  FO 
GO  TO  226 

223  DO  224  1*1,6 

224  SUH  *  SUM*COEFF(K, I ♦ I ) *FO»* I 
PRINT  220,  SUM 

END  OF  ERROR  ANALYSIS 

226  IFISELI-I . 149,47,49 

47  IFl TAU*ZNUM*DTAU <6,46,46 

ACCELERATED  GAUSS-SEIDEL  SOLUTION  USING  A  CONSTANT  ACCELERATION 
FACTOR  AND  WEGSTEIN  ACCELERATION 

61  DO  62  N* I .NODES 
TTtN, I )  *  TIN) 

62  TOLD  IN, I )  *  TIN) 

1  APPLY* I  STOP »NOGS 

63  KSWPS  *  KSWPS* 1 

INDEX  *  XSWPS-KSWPS/3»3 
IFl INDEX  164,64, 65 

64  J  «  3 

GO  TO  66 

65  J  *  INDEX 
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>A 


£ 


a 
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66  IFlKSHPS- ITERS >87,67,28 

67  DO  88  N* I, NODES 
IFlN-l ) 89, 69, 68 

68  IFlN-NODES)72,7I ,7! 

69  FT( I ,3) * (A*TF*Cl I )«TPR( I  )*Dl  I )*TT(2,3) )/t A*Cl I > «Dl l)> 

GO  TO  73 

71  FT I NODES ,  3 )  *  I C I NODES ) *  TPR ( NODES  >  *D I NODES  *  1 1 *TT  < NODES- 1 , 3 ) )  / 1 C I NODE 
IS) *D INODES- I t ) 

GO  TO  73 

72  FT(N,3> • ICIN) *TPR(N>  »0(N- 1 ) *TTlN- 1 ,3) *B(N>  *TT(N* I ,3> )/lCIN}4BlN- l > 
I ♦DIN) > 

73  TEMP*ABSiFTlN*3)-TTlN*3> ) 

IF(TEMP-ZLIM)74,74,75 

74  MOUNT  »  MOUNT ♦ I 

I F ( MOUNT-NODES  >82,37,97 

75  MOUNT  *  0 

82  1F( 3-3)76,79,79 

76  TTIN.3)  »  FTIN.3) 

TTIN,3« l )  »  FTlN,3> 

TOLD1N, 3* ! )  «  FT(N,3> 

IFlKSHPS-2)S3,77,77 

77  IF(KSWPS-I$T0P>78,83,83 

78  TTiN.JM)  *  TT(N, 3)  ♦FACTOR*  ( TOLDlN,  3)  -  TT(N,3) ) 

T0LD(N,3*M  *  TTtN, 3* I ) 

GO  TO  83 

79  TTtN, 3)  •  FT l N , 3 > 

TTtN, I )  •  FT(K,3) 

TOLD ( N , I )  *  FT(N,3) 

IFIKSKPS- ISTOP)8i ,83,83 

81  TT(N.I)  *  TT(N.3>*FACT0R*IT01D(N,3)-TT(N,3> ) 

TOLD5N.S)  *  TTtN, I ) 

83  CONTINUE 

l F ( KSWPS - 1 APPLY  >63,85,85 

85  DO  36  N* I, NODES 
IF(3-2)86,87 ,88 

86  ALPHA  «  1 TOLO(N,2 J -TQLDtN, I > )/( TOLDlN, t ) -TOLD(N,3) > 

GO  TO  89 

87  ALPHA  *  (TOLDlN, 3) -TOLD(N,2 > )/( TOLDIN.2 > -TOLDlN, I > > 

GO  TO  89 

88  ALPHA  »  (TOLDlN, 1 ) -TOLDlN, 3) >/(TOLD<N,3> -TQLDtN, 2) ) 

89  I F I  ALPHA ) 96 , 96 , 9 i 

91  IF  I  ALPHA- 1 , >92,96,96 

92  Q  *  ALPHA/ 1  ALPHA- I . ) 

IFl Q* 1 00 • )99, 1 01 , 1 0 1 

99  Q  *  100. 

101  IF  1 3-3 1 94, 95, 95 

94  TTIN,3‘ I >  *  TT(N,3;*Q»(T0LDIN,3)-TT(N,3) ) 

TOLDlN, 3* I >  «  TTtN, 31 

GO  TO  96 

95  TTIN.I)  •  TT(N,3) *0* I  TOLD! N, 3 1 -TTIN.31 > 

TOLDlN, 1)  •  TTIN.I > 

98  CONTINUE 

{APPLY  *  IAPPLY* INTER 
GO  TO  63 

97  TT(N,3)  »  FT(N,3) 

DO  98  N*l, NODES 
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88  TtN»  »  TT1N.U) 

C  END  OF  ACCELERATED  ITERATIVE  SOLUTION 
GO  TO  24 
C 

C  TRIDAG  SOLUTION 
35  TAU*TAU»OTAU 

SK*I . ‘TAU»20./$2NUH*DTAUI 

K»SK 

RK*K 

I F ( TAU - ZNUK  *DTAU  >251 ,252,46 

251  TF*TSIK1 » ( 7SSK* I ) -TSIlO  >«!$K-RIO 

HF I LM-HFVAR t K ! * ( HFVAR I K* I » -HFVAR ( K I » • ( SK-RK I 
GO  TO  253 

252  TF*TS«K) 

HFILM*HFVAR(K» 

253  PRINT  20,HFILM|TF 

A*  I  ./III  ./HFILMI  ♦DELXm/!24.*ZKC  i  m 
AR( I )*0. 

8RI I  f  *A*Dl  I  I  »CI  l  i 
CRI I )«-D! ! J 

DR! ! >  *Ct  1 1 »TPRI 1 1 *A*TF 
N* NODES 
ARlN>*-D(N-H 
BRlNJ«DlN-ll*CtN» 

CRINJ‘0.0 
DR(N1*C(N)*TPR<NI 
NP* NODES -I 
DO  31  N*2,NP 
ARlN) «-D!N- I ) 

BRiN)»DiN-IJ‘D»N)*C(N) 

CR(N) *-D(NJ 
3!  DRiNl»C(Nl*TPRlH) 

BETA! I ! *BR ! I  I 

GAMMA! I  I *DR( 1 J /BETA!  I  ) 

DO  22  N»  2. NODES 

BETACNJ  *BR!N) -ARINI »CR!N* I J /DETAIN- 1 1 

22  GAMMA! N) « iDR IN) -AR IN) ‘GAMMA! N* I > >/BE?A(N> 

N* NODES 

TNEWIN) ‘GAMMA! Nl 
LAST ‘NODES- 1 
DO  23  L* I .LAST 
N'NODES-L 

23  TNEKt  N) « GAMMA !NJ  -CRIN5  *TNEMlN*  t  I/8ETAIN) 

44  N*NODES 
PRINT  270, TAU 
PRINT  240 

PRINTSO,  <N, TNEWIN I ,N* I .NODES > 

KOUNT*0 

32  DO  34  N*1 .NODES 
34  TPRINJ  «TNEN(N) 

I F  t  NODES  I -NODES  1 43 , 1 38 , 46 
49  IF! TAU-ZNUM»DTAU>35 ,46.46 

45  CALL  EXIT 
END 
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LISTING  OF  THE  COEFFICIENTS  FOR  THE  CHEBYSHEV  POWER  SERIES 


COEFFICIENTS  FOR  SPATIAL  TRUNCATION  ERROR  CURVE,  8IOT  NO,*O.S 


10085509E»02  ,3t3!3865E»02 
«6G5768S5£*0 i  I429433IE»0l 
•.30IS2627E-94 


-,40438398E*02  ,32S76778E»02  -. 17I9S332E*02 

.22260 I58E*00  -.21906I65E-0I  . I2322493E-02 


COEFFICIENTS  FOR  SPATIAL  TRUNCATION  ERROR  CURVE,  BIOT  NO.«l.O 


-.4774I800E»0!  . I07I6588E*02 

•.636968201*01  , 1 8?48795E*0 1 

.707S2864E-04 


,39482730E»0!  - .  I S347922E-Q2  . !337898S£*02 

• ,349031 0 IE*00  .3997IS59E-0I  • .2S676872E-02 


COEFFICIENTS  FOR  SPATIAL  TRUNCATION  ERROR  CURVE,  BIOT  NO. *3.0 


-.! 074266 l £*0?  .S79S30S2E*02 

. I2I3S378£*02  « .29680640E*0 I 
-.7373634IE-04 


.78I97286E-02  .628073S2E*02  - .33S33087£*02 

•48I0!OI7E*00  - .49327886E-0 l  .289I9389E-G2 


COEFFICIENTS  FOR  SPATIAL  TRUNCATION  ERROR  CURVE,  BIOT  NO. *10.0 


*.473436S0E*0I  .57069I20E*02 

. I585902IE»02  -.41 I68I94E*0I 
-. I20440I4E-03 


.84933970E*02  .7242384QE*02  -.41 I2I90SE»02 

.70284882E*00  -.7S32050IE-0I  . 458 1 0 1 70E- 02 


COEFFICIENTS  FOR  SPATIAL  TRUNCATION  ERROR  CURVE,  BIOT  NO, •INFINITY 


*. !39050I2£»02  • I 065I562£*03 

,35326802E*02  -,83I97870E*0I 


-.174961 I7E03  . I5979799E*03  • ,92659020E*02 

. I4726003E*0I  - • I 5254262E*00  .896656I3E-02 


• .22953966E-03 


C  COEFFICIENTS  FOR  TIME  TRUNCATION  ERROR  CURVE,  BIOT  NO. *0.5 

I4942700E»00  . I3063356E«0 I  ,7I703S53E*00  - ,27S29I48E*00  .397I8724E-0I 

-.26482984E-02  .6752I740E-04 


C  COEFFICIENTS  FOR  TIME  TRUNCATION  ERROR  CURVE,  8IOT  NO. >1.0 

-.52064I00E-0I  ,2I0309I5E*0I  .94S643I0E»00  * .4 1 0SI 7!6E»00  .63S6833IE-0! 

-.45599I69E-02  . I2669898E-03 

C  COEFFICIENTS  FOR  TIME  TRUNCATION  ERROR  CURVE,  BIOT  NO. *3.0 

-.674I5600E-0I  .37673700E-0 1  . I466234IE*0I  - .87I03I5IE-00  . I6056679E-00 

I3017263E-01  .3SS4I555E-03 


C  COEFFICIENTS  FOR  TIME  TRUNCATION  ERROR  CURVE,  BIOT  NO. *10.0 

-.978I5400E-0I  ,S05!0I20E*0I  . I !366794E*0 I  -.8608I270E-00  . I6645989E»00 

-. 13872656E-0I  .43I4869IE-03 


C  COEFFICIENTS  FOR  TIKE  TRUNCATION  ERROR  CURVE,  BIOT  NO. *  INFINITY 

16628430E»00  .647592 1 0E*0 1  .49463530E-00  - ,7488!550£»00  . IS6669I0E»00 

-. I3388434E-0I  .4I87783IE-03 

These  coefficients  are  used  in  storing  the  spatial-  and  the  time- 
truncation  error  curves  in  the  ONE-D  program. 
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Section  2 


ANALYTICAL  SOLUTION  FOR  AN  INFINITE  PLATE 


C  SOLUTION  OF  TEMPERATURE  DISTRIBUTION  IN  A  FINtTE  SLAB  ~  ONE  DIMENSIONAL 
C  FLOW  -  FINITE  FILM  COEFFICIENT  ON  ONE  FACE.  THE  OTHER  FACE  INSULATED 
C  ALPHA  *  THERMAL  DIFFUSIVITY,  SQ.FT. /HR. 

C  AM  «  EIGEN  VALUES  OF  M*TAN!M) 'CONSTANT 
C  BIOT  •  BIOT  NUMBER 

C  DELS  ■  LENGTH  OF  PLATE  FROM  SURFACE  TO  INSULATED  FACE.  INCHES. 

C  DELX  •  SPATIAL  INCREMENT,  INCHES 
C  DTAU  »  TIME  STEP,  SECONDS 
C  HFILM  •  FILM  COEFFICIENT,  BTU/HR-SQ.FT, -F, 

C  TAU  •  TIME,  SECONDS 
C  TF  •  FILM  TEMPERATURE,  F. 

C  Tt  *  INITIAL  TEMPERATURE,  F. 

C  THETA  «  FOURIER  NUMBER 

C  X  »  DISTANCE  FROM  INSULATED  FACE  TO  NODE  CONSIDERED,  INCHES 
C  ZC  «  SPECIFIC  HEAT.  6TU/LBH-F. 

C  ZK  *  THERMAL  CONDUCTIVITY,  BTU/HR. -FT. -F. 

C  ZRHO  «  DENSITY.  LBH/CU.FT. 

C 

DIMENSION  AMI 59 ) .OEXPISO I .SINS  150 ) ,SIN2!50 ) 

10  FORMAT I7FI0.0.2S5I 

20  FORMAT 1 1  HO , 1  OX , 43HTRANS I ENT  SOLUTION  FOR  INSULATED  FLAT  PLATE) 

30  FORMAT! IHO , 1 0X.2HK* ,F7.3, I0X.2HC*  »F?.4,5X,8HDENSITY*,F9.3) 

AO  FORMAT! IH  .I0X.I8HFILM  COEFFICIENT* ,FI0. 31 
50  FORMAT! IK  , I 0X.6HALPHA* ,F8. 4,2X,2HL* ,P8.4 > 

60  FORMAT! IH  , ) OX. I3HINITIAL  TEMP* ,F 1 0 . 4,2X, I0HFILH  TEHP*,FI0.4I 
70  FORMAT !8F 1 0.0) 

80  FORMAT  1 1  HO . 1 2X . 1 2HE I  GEN  VALUES) 

SO  FORMAT! IH  , 10X.2HH! , I3.3H) •  .FI0.5) 

100  FORMAT! FI 0.0) 

110  FORMAT 1 1  HO . 1 2X , 4HT I  ME ■ 8X , 3HX/L , 4X , I OHTEMP  RATIO.GX.4HTEKP) 

120  FORMAT! IH  , I 0X,F3.2,3X,F7.4,4X,FS.6,5X,FS.3) 

130  FORMAT 12H0, 1 0X.5HDELX* ,F7.4, 3X, I0HBIOT  NO.  » .FS.4.3X, I2HFOURIER  NO 
I . * ,F8.4) 

READ  10. ZK.ZC. ZRHO, HFILM, DELS. TF.TI.KODE.N 

PRINT  20  \ 

PRINT  30. ZK.ZC, ZRHO  i 

PRINT  40, HFILM  ) 

BIOT*HFILH*DELS/l I2.9*ZK) 

ALPHA «ZK/!ZC* ZRHO)  ! 

PRINT  50, ALPHA, DELS  | 

PRINT  60, TI ,TF  j 

PRINT  80  i 

IFIKODE- I  I  2,1,2 

1  READ  70,1  AM <I),I-*t,N) 

GO  TO  3 

2  CALL  EIGEN(N,AM,8IOT) 

3  DO  4  1*1, N 

4  PRINT  90, 1, AMID 
READ  1 00, DELX 

IF! SENSE  SWITCH  I)  IS, 5 

5  READ  1 00. TAU 
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THETA*ALPHA*TAU/ (3800.0 ‘DELS ‘DELS/ 1 44 . > 

PRINT  130, DELX, BIOT. THETA 
00  8  I«!.N 
SINK  l )«SINFIAH(  I  i I 
S1N2( I )*SINF(2«0*AM( I 1 I 

6  0EXP( M *£XPF( «AH( I ) **2*TH£TA) 

X*0, 

7  SUM*0. 

00  8  I*1,N 

8  SUM*SUMMSIN1 1  n*DEXPl !  »/C2,0*AH(  !l*SIN2t  I »  M*CO$F(AK(  II*X/DEL$1 
TRATIO*4iO*$UH 

TEHP*TF-»(TI-TFI*TRATIQ 

XRATIQ-X/DELS 

PRtNT  110 

PRINT  I20.TAU.XRATIO.TRATIO.TENP 
X«X»DELX 
IF(X*DELSi  7,7,5 
15  CALL  EXIT 
END 
C 
C 

SUBROUTINE  EIGEN(N.AH,8lOTf 

DOUBLE  PRECISION  T,Tl ,ANeZK,2C,2RHO,TF,HF!LK,X,DELX.DELS,A,TAU,DE. 
18107 

C  HALF  INTERVAL  SEARCH  FOR  ROOTS  OF  COT(H! -H/BIOT-O. 

DIMENSION  AM(50) 

PI  *3.  Ml 59265 
EPS* l . CE-3 
1*0 

1  I  - 1  ♦  I 

I EYE  *1-1 

EYE*IEYE 

A«?l/I80.*EYE*Pl 

8*Pl/2. ♦EYE*PI 

FA«DCOS( A J/DSINIAI -A/BIOT 

2  X*  t A»8 1/2. 

COT  *DCOS ( X I /OS  I N ( X I 
XOB*X/BlOT 
F-COT-XOB 
I F ( F I  3,11,4 

3  IF(F*EPSJ  5,11,11 

4  IF(F-EPS)  11,11,5 

5  IF(F»FAI  6,11,7 

6  B*X 
FB«F 

60  TO  2 

7  A*X 
FA*F 
GO  TO  2 

It  AH( I ) »X 

!F(N-l>  12,12,1 
12  RETURN 
STOP 
END 
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MODIFIED  HEAT-TRANSFER  PROGRAM  "ONE-D”  WITH 
ERROR  ANALYSIS  AS  APPLIED  TO  RAMP-FUNCTION 
BOUNDARY  CONDITIONS 

The  program  used  to  generate  data  for  the  ramp-function  boundary 
condition  is  an  abbreviated  form  of  ONE-D.  Since  the  TRIDAG  method  of 
solution  was  the  fastest  and  simplest  to  use,  it  was  chosen  and 
revised  to  handle  only  homogeneous  solids.  It  was  also  revised  so 
that  the  automatic  time  step  generator  subroutine  (ATSG)  could  be 
used. 

The  error  routine  portion  of  ONE-D  is  also  incorporated.  Revisions 
were  made  in  this  routine  so  that  the  spatial  error  for  either  time 
step  functions  or  ramp  functions  is  given.  Time-truncation  error  is 
printed  only  for  step-function  boundary  conditions. 


CONTROL  VARIABLES 

K0DE(1)  If  K0DE(1)  is  1,  the  automatic  time-step  generator  is  used. 
Values  for  ANUMBE,  PRCNTH,  and  PRCNTL  must  be  read. 

PRCNTH  and  PRCNTL  are  given  in  decimal  form, 

KODE(2)  If  K0DE(2)  is  1,  a  Fourier  number  is  read  and  the  time 

step,  DTAU,  is  calculated.  If  K0DE(2)  is  not  1,  DTAU  is 
read, 

K0DE(3)  If  KQDE(3)  is  1,  the  initial  temperature  of  the  body  is 

set  at  zero  and  the  time,  TAU,  is  set  at  zero.  If  K0DE(3) 
is  not  1,  values  for  the  temperature  at  the  nodes  and  TAU 
are  read. 

KODE(4)  If  K0DE(4)  is  1,  the  temperature  of  the  nodes  may  be  found 
for  a  number  of  specified  times.  If  K0DE(4)  is  not  1,  the 
temperatures  will  be  printed  according  to  DTAU. 

TAU1  TAU1  may  be  used  as  a  control  variable  in  addition  to 

specifying  the  time  at  the  end  of  the  ramp.  If  TAU1  is 
less  than  1,  a  step  function  will  be  assumed. 

TAU2  may  also  be  used  as  a  control  variable.  TAU2  is 
ordinarily  used  to  specify  the  time  at  which  the  run  will 
stop.  However,  if  TAU2  is  0,  then  the  user  can  specify 
the  terminal  temperature  in  any  node  desired. 


TAU  2 


ARRANGEMENT  OF  INPUT 

Data  Card  1; 

FORMAT  (715),  KODE(l) ,  KODE(2) ,  KODE(3)  ,  KODE(4) ,  NODES, 

LTP,  NTP.  This  card  must  be  read  for  each  set  of  data. 

Data  Card  2: 

FORMAT  (3F10.3),  TAU1,  TAU2,  FT.  This  card  must  be  read 
for  each  set  of  data. 

Data  Card  3: 

FORMAT  (F1Q.3,  15),  TSTOP,  NSTOP.  This  card  is  read  if 
TAU2  is  0, 

Data  Card  4: 

FORMAT  (3F10.3),  ANUMBE,  PRCNTH ,  PRCNTL.  This  card  is  read 
if  K0DE(1)  is  1. 

Data  Card  5: 

FORMAT  (4F10.3,  E10.3,  15),  DELX1,  ZK1,  ZRHQl,  ZC1,  ZH1,  M2. 

This  card  must  be  read  for  each  set  of  data. 

Data  Card  6: 

FORMAT  (E10.3)  HFILM.  This  card  must  be  read  for  each  set 
of  data. 

Data  Card  7; 

FORMAT  (F10.3)  FO.  This  card  is  read  if  K0DE(2)  is  1. 

Data  Card  8: 

FORMAT  (F10.3) ,  DTAU .  This  card  is  read  if  K0DE(2)  is  not  1. 
Data  Card  9: 

FORMAT  (2F10.3)  TAU,  TF.  This  card  is  read  if  K0DE(3)  is  not  1. 
Data  Card  10-22: 

FORMAT  (8F10.3),  TPR.  These  cards  are  read  if  KODE(3)  is  not  1. 
The  number  of  data  cards  is  dependent  on  the  number  of  nodes  in 
the  system  with  a  maximum  of  100  nodes. 
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Data  Card  23: 


FORMAT  (12) ,  NO.  This  card  is  read  if  KODE(4)  is  1. 
Data  Card  24-25: 


FORMAT  (8F10.3).  These  cards  are  read  if  KODE(4)  is  1. 

Data  Cards  26-56: 

FORMAT  (5E15.8).  These  cards  are  used  to  supply  the  coefficients 
of  the  power  series  that  describe  the  spatial-  and  time-truncation 
error  curves.  If  more  than  one  set  of  data  is  read  (LTP 
greater  than  1)  these  cards  are  read  only  once. 


PROGRAM  LISTING 


C  TRIDAG  SOLUTION  OF  ONE-DIMENSIONAL  TRANSIENT  HEAT  TRANSFER 
C  PROBLEM  KITH  ERROR  ANALYSIS 
C 

C  NOTATION 

C  ANUKBE  *  THE  NUMBER  OF  TIME  STEPS  BETWEEN  CHECKS  BY  ATSG 
C  AR.BR.CR.DR  *  COEFFICIENT  ARRAYS  CONTAINING  THE  SUB-DIAGONAL, 

C  DIAGONAL.  SUPER -DIAGONAL,  AND  RIGHT  HAND  ELEMENTS 

C  OF  THE  TRIDIAGONAL  SYSTEM 

C  BIOT  3  BIOT  NUHBER 
C  DELX  3  SPATIAL  INCREMENT,  INCHES 
C  DTAU  •  TIME  STEP,  SECONDS 
C  FO  «  FOURIER  NUMBER 

C  FOI  «  FOURIER  NUMBER  AT  THE  END  OF  THE  RAMP 
C  FT  «  FINAL  FLUID  TEMPERATURE,  F. 

C  HFILH  *  FILM  COEFFICIENT.  BTU/HR. -SQ.FT. -F. 

C  KODEINI  *  CONTROL  VARIABLES 

C  LTP  *  NUMBER  OF  SETS  OF  DATA  TO  BE  READ 

C  M23  NUHBER  OF  NODES  BETWEEN  ZERO  THICKNESS  NODES 

C  NO  *  NUMBER  OF  TIMES  AT  WHICH  A  TEMPERATURE  PRINT-OUT  IS  REQUESTED 

C  NODES  *  TOTAL  NUMBER  OF  NODES  IN  THE  SYSTEH 

C  NTP  3  NUMBER  OF  TIME  STEPS  BETWEEN  TEMPERATURE  PRINTINGS 

C  NSTOP  *  NODE  SELECTED  TO  TERMINATE  RUN 

C  PRCNTH  *  LARGEST  PERCENT  INCREASE  IN  TEMPERATURE,  EXPRESSES  AS  k  DECIMAL 
C  PRCNTL  *  SMALLEST  PERCENT  INCREASE  IN  TEMPERATURE.  EXPRESSED  AS  A  DECIMAL 
C  SCOEFF  3  COEFFICIENTS  OF  POWER  SERIES  FOR  SPATIAL  ERROR  ANALYSIS 
C  TAU  3  TIME,  SECONDS 

C  TAU 1 3  TIME  AT  THE  END  OF  THE  RAMP,  SECONDS 
C  TAU2  3  TIME  AT  TERMINATION,  SECONDS 
C  TAUT  3  TEMPORARY  TIME  STORAGE 

C  TCOEFF  3  COEFFICIENTS  OF  POWER  SERIES  FOR  TIME  ERROR  ANALYSIS 
C  TF  3  FLUID  TEMPERATURE  AT  SPECIFIC  TIME  UNDER  CONSIDERATION,  F. 

C  TIMEUI  3  TIMES  AT  WHICH  A  TEMPERATURE  PR  I  NT -OUT  IS  REQUESTED,  SEC. 

C  TNEW  •  TEMPERATURE  OF  NODE  AT  END  OF  ITERATION  OR  TIME  STEP,  F. 

C  TPR  3  TEMPERATURE  OF  NODE  AT  PREVIOUS  TtME  STEP.  F. 
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C  TPRT  *  TEMPORARY  STORAGE  FOR  PREVIOUS  TEMPERATURES 
C  TSTOP  *  TERMINATION  TEMPERATURE  FOR  NODE  NSTOP,  F. 

C  ZC  *  SPECIFIC  HEAT,  BTU/LBH-F. 

C  ZH  »  CONDUCTANCE  8ETWEEN  NODES.  TO  ALLOW  FOR  CONtACT  RESISTANCE 
C  ZK  «  THERMAL  CONDUCTIVITY,  8TU/HR. -FT. *F. 

C  ZRKO  *  DENSITY,  L8M/CU.FT. 

C 

DIMENSION  KODEtSJ ,TPR< (GO  ) ,TNEN( 1001 ,C( 1001 ,CRC 100 ( .0(100) ,AR( 100) 
l ,8R( 100) ,DR( 100) ,BETA( 1 00 1 , GAMMA) 100) ,ZK( (00 ) ,ZC( 100) , ZRKO 1 100 ) , 
2ZH( 1 00 ) ,OELX( 100), TPRT ( 1 00) , TINE! 10! , SCOEFF 1-7,111, TCOEFF ( S , 7 1 
COMMON  TNEK , TPR , TAU ,DTAU , NODES , TF . FT , TAU l . NUMB . TPRT . DTAUT . TAUT . 

1 PRCNTH .PRCNTL ,C,ZC, ZRKO , DELX 
C 

10  FORMAT ( 7 1 S ) 

20  FORMAT (3F l 0,8) 

30  FORMAT (8F 1 0.3) 

40  FORMATiFIO.3) 

SO  FORMAT ( IH) , 4HNODE , I IX.7HDELTA  X. MX, IHK, ISX.7HDENSITY. ISX, 1HC, ISX, 
12HHC) 

60  FORMAT(T2, 13,4(5X,FH.3) ,9X,£I0.3> 

70  FORMAT! IH0.8HDELTA  T*,F12.4,SH  SEC.) 

80  FORMAT! Tt ,5F! i ,3,5X,5Ft I ,3) 

90  FORMAT! EJ 0.3) 

1 00  FORMAT  <1H0, I3HHF l LM*CONST . « , F 1 2 . 3 ) 

HO  FORMAT i ) HO ,  1 2HFLUID  TEMP. «, FI  0.3) 

120  FORMAT (12) 

130  FORMAT! IHO, 1 1X.4HTIME.F12. 4, 1  OK  SECONDS) 

HO  FORHAT(4Fi0.3.EI0.3.IS) 

150  FORHAT12FI0.3) 

ISO  FORMAT! IHO,  .MAXIMUM  SPATIAL  TRUNCATION  ERROR  IS  LESS  THAN  -2.0 
1  PERCENT.) 

170  FORMATI5EI5.8) 

180  FORMATIFIO.3, IS) 

190  FORMAT! 1H0H9HMAXIMUM  SPATIAL  TRUNCATION  ERROR  IS  *2.5  PERCENT  ) 
200  FORMAT! IH0.37HMAXIMUM  SPATIAL  TRUNCATION  ERROR  IS  -.F6.1.9H  PERCEN 
IT  ) 

210  FORMAT! I HO. 77HT I ME  STEP  IS  TOO  BIG.  SUGGEST  USING  SMALLER  FOURIER 
INUMBER.  FOURIER  NUMBER  IS.  FI0.2) 

220  FORMAT (IHO, 34HMAX l MUM  TIME  TRUNCATION  ERROR  IS  -.FS.I.9H  PERCENT  ) 
230  FORMAT! 1H0,4SHMAXIMUM  SPATIAL  TRUNCATION  ERROR  IS  *1.3  PERCENT  ) 

C 

JKT«0 

KKK*0 

I  READ!  105,10)  NODE!  I  )  ,KODE(2)  ,KODE(3)  .KODEH)  .NODES  LPT.NTP 
READ (5, 20)  TAUI.TAU2.FT 

IFI TAU2.E0. 0 . )  READU05.I80)  TSTOP, NSTOP 
I F ( MODE ! 1 ) « EQ . 1 )  READ! 105,20)  ANUMBER, PRCNTH, PRCNTL 
NUMBER* ANUM8ER 
NUMB*0 

READ! 105, HO)  DELX1 ,AKI , ZRHOl , ZC 1 ,ZH1 ,M2 
DO  42  N* 1, NODES 
DELX(N) *DELXI 
ZKlNHZK! 

ZRHOCN) *ZRKOi 
ZHIN)*ZH1 
42  ZC(N) *  ZC 1 

DO  3  N* I .NODES, M2 
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3  DELXiN) *0. 

WRITE! 108^50) 

DO  4  N*l. NODES 

4  WRITE! 108,60)  N.DELXtN) ,ZK(N> ,ZRHO(N) ,ZC(N) ,ZH!N) 

READ! 105,90)  HFILH 

WRITE*  108, 100)  HFIL.N 
IF!KODE(2)*l >  8,7,8 

7  READ! 105,40)  FO 

DTAU*FO»DELXI •«2*2RHOI ‘2CI •6.25/ZKI 
GO  TO  9 

8  READ! 105,40)  DTAU 

9  WRITE! 108,70)  OTAU 
DTAUT«DTAU 
IF(TAUI)  19,13,18 

13  TF»FT 
GO  TO  17 

16  IF(DTAU.GE.TAUl)  GO  TO  13 
7F*FT«DTAU/TAUI 

17  WRITE!  108,.!  10)  TF 

! F ! KOOE 1 3 ) - 1 >  14,12,14 

14  READ! 105,150)  TAU,TF 
WRITE! 108, 130)  TAU 

READ! 105,30)  ! TPfttN) ,N« ! .NODES) 

WRITE! 108,80)  (TPR(N) ,N* I , NODES) 

GO  TO  18 

12  DO  IS  N* I, NODES 

15  TPRlN) *0, 

TAU»Q . 

18  KOUNT*0 
NP* NODES -l 
DO  19  N* 1 ,NP 

19  D(N) *  I ./ (DELXiN)/ (24. *ZK!N) ) ‘DELXIN* I )/!24.  »ZIC(N* 1 ) I ♦ I »/ZH(N) ) 
DO  39  N* I .NODES 

TPRTlN) *TPRlN) 

39  C ! N ) * 30 0 . • ZC ! N ) »ZRHO I N I *DELXlN ) /DTAU 
A* 1 • / ( 1 ,/HFILH) ♦DELX! I )/!24.*ZK< I ) ) ) 

I F ! KODE ! 4 ) *  1 )  35,95,35 
95  READ! 105, 120)  HO 

READ! 105,30)  i TIME (N) ,N* I ,NO) 

TAUT * TAU 
NN*  I 

35  NUMB * NUMB ♦! 

TAU«TAU*DTAU 
I F ( KODE ! 4 ) ' I)  83,61,83 

61  IF! NN.EQ.NO* 1 )  GO  TO  63 
IF! TAU-TIMElNN) )  63,62,62 

62  TAU* TIME  INN) 

NN*NN«! 

IF! TAU-TAUI )  64,65,65 

64  TF*FT*TAU/TAUI 
WRITE! 108, 1 10)  TF 
GO  TO  63 

65  TF*FT 


i 
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63  MOUNT *KOUNT • 1 
AR( I ) *0. 

BRm»A»Dtl)*CU> 

CR( 1 )»-D( I ) 

QRC  n*CU>*TPRm*A*TF 
N; NODES 
AR(N) *-DtN* I ) 

BR(N) »D(N* I ) *C(N) 

CRtNl *0 . 

DR!N)«C!N>*TPR(N) 

DO  21  N*2,NP 
AR!N)*-D(N-l ) 

8R(N) *0!N* I ) *D(NJ  *CJN) 

CR(N) *-D!N) 

21  DR(N) *C!N) *TPR(N) 

BETA(N*BR!  1) 

GAMMA! 1 )»DR( I > /BETA! I ) 

DO  22  N*2, NODES 

8ETA!N)«SRIN)*AR(N)*CR(N-1)/0ETA<N-I ) 

22  GAMMA!N)*(DR(N>~AR(N)*GAMMA(N*i > 1/BETAlN) 

TNEWt NODES) ■GAMMA! NODES) 

DO  23  L* I ,NP 
N*NODES-L 

23  TNEWlN) * GAMMA ! N ) -CR ( N ) • TNEW ( N ♦ I ) /BETA ( N ) 

IF! MOUNT. LT.NTP)  GO  TO  33 

28  WRITE ( E 08,130 )  TAU 

MR  1 IE ll 08 , 80 )  (TNEWlN) ,N«1 , NODES ) 

MOUNT* 0 

C  SPATIAL  AND  TIME  TRUNCATION  ERROR  ANALYSIS 
IF(MMM.GE.I)  GO  TO  229 

C  READ  COEFFICIENTS  OF  POWER  SERIES  FOR  SPATIAL  TRUNCATION  ERROR 
DO  .*89  M*  1 , 7 

199  REAL* I  OS, 1701  « SCOEFFlM, I ) . I  *  I , II 
229  BIOT*HFILM»bELXI/24./ZKl 

FO«ZMI *TAU/(6.25*ZC1 «ZRNOI «0ELXl ••2) 

C  DETERMINE  IF  RAMP  FUNCTION  OR  STEP  FUNCTION 
IFlTAUi.LT.l.)  GO  TO  222 
FOI »ZM! «TAU1/(6.25«ZCI »ZRH01 •DELXI *«2) 

IF(FO.LE.FOl)  M*6 
IF(FO.GE.FOl)  K*7 
IFIFO.LE..7S)  M*6 
IF(FO.GT. 10. )  M*6 
IFtFO.LE.20. )  GO  TO  226 
WRITE! 108, 160) 

GO  TO  225 

226  IFtFO.GT. .25)  GO  TO  217 
WRITE! 108. 190 ) 

GO  TO  225 

C  DETERMINE  SPECIAL  CONDITIONS 
222  IFtrO.LE.7. >  GO  TO  202 
WRITE! 108.230) 

GO  TO  225 

202  IFIBIOT.GT. .5)  GO  TO  205 
IFtFO.GE. .75)  GO  TO  216 
WRITE! 1 08. 190 ) 

GO  TO  225 
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205  IFlFO.GE. .25)  GO  TO  207 
WRITES  108.190) 

GO  TO  225 

C  DETERMINE  WHICH  CURVE  TO  USE,  BASED  ON  BIOT  NUMBER 

207  IFIBIOT-IO.)  211,209,208 

208  K  =  5 

GO  TO  217 

209  K  =  4 

GO  TO  217 

21  I  IFlBlOT-3. )  213,212,209 

212  K*3 

GO  TO  217 

213  IF! BIOT- I . )  215,214,212 

214  K *2 

GO  TO  217 

215  IFtBIOT .GT . .5)  GO  TO  214 

216  K*  I 

217  IF(TAU.GT.DTAU)  GO  TO  218 
IFII.LE.II)  60  TO  223 

218  SUM  =  SC<j£FF ( K,  I  ) 

DO  2 1 9  1*1.10 

219  SUM=SUM-SCOEFFt K, I  * ! ) «FO«» I 
I F l SUMLT . 0 .  !  SUM* -SUM 
WRITE! 108.200)  SUM 

225  IFtKKK.GE.  I  )  GO  TO  228 
KKK=KKX‘I 

C  READ  COEFFICIENTS  OF  POWER  SERIES  FOR  TIME  TRUNCATION  ERROR 
DO  221  K*l  ,5 

221  READ! I0S, 170)  ! TCOEFFIK, I ) , I  *  1 ,7 > 

228  IF! TAUI . GE. I . )  GO  TO  33 
IF(TAU.GT.DTAU)  GO  TO  33 
FO*ZKI •DTAU/(6.25*ZC1 •ZRHOl»DELXi ••2) 

I  FIFO. LE . 10. )  GO  TO  207 
WRITE! 108,210)  FO 
GO  TO  33 

223  SUM'TCOEFF ! K, I ) 

DO  224  1=1 ,6 

224  SUM*SUM»TCOEFF (K, I  *  I ) »FO* • J 
WRITE! 108.220)  SUM 

33  IFIKODE! 1 ) -  I  I  38,32,38 
32  IF! NUMB -NUMBER  1  38,37,37 

37  CALL  ATSG 
GO  TO  29 

38  DO  24  M* I, NODES 

24  TPR t N ) «TNEW(Ni 

I F ( TAU-TAU I )  25,29,29 

25  TF  =FT» ( TAU*DTAU I /TAUI 
IF(TF.GE.FT)  TF*FT 

29  IFtKODE! 41  - 1 )  67,66,67 

66  IF ! TAU-T I  ME ( N/ ) )  31 ,67,67 

67  IF ( TAU2 . EQ . 0 . )  GO  TO  68 
IF ( TAU. GE . TAU2 )  GO  TO  26 

31  WRITE! 108, I  10)  TF 
GO  TO  35 

68  I F I TNEW ( NSTOP ) . GE. TSTOP I  GO  TO  26 
GO  TO  31 
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26  JKT  « JKT ♦ l 

IFlJKT-LTP)  1,41.41 
41  STOP 
END 


AUTOMATIC  TIME-STEP  GENERATOR  SUBROUTINE 


SUBROUTINE  ATSG 

DIMENSION  TPRI I00I.TPRTI  100 ) .TNEW! 100 ) ,C( 1 00 ) ,ZCI 1 00 > .ZRHOl I QO ) , 
IDELXUOO) 

COMMON  TNEW , TPR , TAU , DTAU , NODES  ,  TF , FT , TAU ! , NUMB , TPRT . OTAUT , TAUT , 

1 PRCNTH . PRCNTL , C , ZC , ZRHO.DELX 

10  FORMAT tl HO, SX.5SHDELTA  TAU  WAS  HALVED.  TIME  SET  BACK  TO  LAST  CHECK 
••POINT.,  1H  ,5X,37HTPR,S  ARE  LISTED  BELOW  FOR  THAT  TIME.! 

20  FORMAT! I  HO. I  OX. 7HT1 ME  IS.F12.S.SX. I6HNEW  DELTA  TAU  IS.F10,5) 

30  FORMAT! 1H0 .5X.26HPREVIOUS  NODE  TEMPERATURES) 

40  FORMAT! 1H  , SF 1 1 .3.SX.5F11 .3) 

SO  FORMAT! IHO.SX.S2HOELTA  TAU  WAS  UNCHANGED.  PROGRAM  PROCEEDS  AS  BEFO 
IRE. ) 

60  FORMAT! I  HO, SX.48HDELTA  TAU  WAS  DOUBLED.  TEMPERATURES  ABOVE  ARE  OK) 
NUMB*0 

MODES*NODES/2*2 
DTUP*PRCNTH«TF 
DTLOW*PRCNTL*TF 
DO  I  N*2, NODES. MODES 
D IFF *ABS! TNEW! N) -TPRlN) I 
IFlDIFF.GE.DTUP)  GO  TO  2 
IFlDlFF.GE.DTLOW)  GO  TO  3 

1  CONTIMUE 
DTAU*DTAU*2. 

WRITE! 108,60) 

WRITE! 108.20)  TAU, DTAU 
GO  TO  4 

3  WRITE! 108,50) 

4  DO  5  N* I, NODES 

C (N) *300 . *ZC( N) *ZRHO(N) *DELXlN) /DTAU 
TPRlN) *TNEWIN) 

5  TPRT(N) *TNEWtN) 

TAUT* TAU 

IFlTAU-TAUI )  6.7.7 

6  TF*FT*!TAU«OTAU)/TAUI 
IF(TF.GE.FT)  TF^F" 

7  RETURN 

2  DTAU*DTAU/2. 

TAU*TAUT 

DO  8  N* I .NODES 

C!N)*300.*ZC!N)*ZRHO!N)*DELX!N)/DTAU 
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8  TPRlN) *TPRT!N) 

WRITE! 108,20)  TAU.OTAU 
WRITE) 108,10) 

WRITE! 108,30) 

WRITE! 108,40)  «TPRlN),N«l .NODES) 
IF<TAU-TAUI)  9,11,11 
8  TF*FT*lTAU»DTAU)/TAUl 
IF! TF.GT.FT )  TF*FT 
RETURN 
II  TF»FT 
RETURN 
END 
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Section  4 

ANALYTICAL  SOLUTION  FOR  AN  INSULATED  FLAT  PIATE 
WITH  RAMP-FUNCTION  BOUNDARY  CONDITIONS 


7  X*0£LX*X 

8  IF(TAU.LE.TAU1 I  CO  TO  18 
3  X*X1 

PROD* ALPHA- ( TAU-TAU I ) /360  0 . / < DELS* -2  I  1 1 44 , ! 

DO  12  U* I .NODES 
$UM*0. 

DO  II  t « 1  ,N 
AS  IN! *DS!N( AM!  I ) ) 

ASIN2*DSIN(2.*AH( t ) ) 

POKER -AM ( U**2*PR0D 
IFIPOKER.GT. 174. I  POWER- 174. 

1 1  SUM* SUM* (AS  IN  I • (DEXP! -POWERI - 1 )/CAMt I ) **2* ( AMI  1 1 *2. *ASIN2I i • 

IDCOSt  AHt  H*X/DELS) ) 

Tl ( Jl *A* ( ( TAU-TAU 1 ) *B«SUM) 

T(JI«T(U)-T1 (J) 

IF(Tl*J) -0.0 1 »  13,13,12 

12  X*X*DELX 

13  TFRAMP*TF/TAUI*TAU 
IFtTF.LT.TFRAMPJ  TFRAMP-TF 

15  WRITEt 108.100)  TAU.TFRAMP 
WRITE! 108, 1 10)  (TOO  ,K* I ,  J) 

TAU*TAU*DELTAU 
IF(TAU-TAU2)  S.5,16 

16  STOP 
END 

C 

C 

SUBROUTINE  EtGEN(N.AM.BlOT) 

DOUBLE  PRECISION  T.TI ,AM,ZK,ZC,ZRHO,TF,HFILM,X,DELX,DELS,A,TAU,DE, 
IBIOT 

C  HALF  INTERVAL  SEARCH  FOR  ROOTS  OF  COT(H»-M/BIOT*0. 

DIMENSION  AM(S0) 

PI  *3.  HI 59265 
EPS* I . 0E-3 
1*0 

1  1*1*1 

I  EYE* l -  I 
EYE* I  EYE 
A*P1/I80.*EYE*PI 
B*Pl/2. *EYE*Pl 
FA*OCOS( A)/DSIN(A) -A/BIOT 

2  X*(A*B)/2. 

COT*DCOS(X)/DSIN(X) 

XOB*X/BIOT 

F-COT-XOB 

IF(F)  3,11,4 

3  IF(F*EPSJ  5,11,11 

4  IFtF-EPSI  11,11,5 

5  IF(F*FAJ  6,11,7 

6  B*X 
F8*F 

GO  TO  2  65 
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C 

C 

c 


1 1 


A*X 
FA*F 
GO  TO  2 
AK!l)*X 


D1UP9 

INSULATED  FLAT 


PLATE,  RAMP  FUNCTION  ON 


TRANSIENT  SOLUTION  FOR 
SOUNDING  FLUID. 

DOUBLE  PRECISION  T.T! , AH, ZK.ZC.ZRHO.TF, HFILK.X.DELX, DELS, TAU, DELTA 
IU.TAUI .TAU2.8IOT, ALPHA, XI .A.B.ASINi , AS IN2, POWER, SUM, TFRAHP, PROD 
DIMENSION  T!200 ) ,T1 (200 ) ,AH!50 ) 

10  FORMAT (AF 1 0.0) 

20  FORMAT ( I K I , I  OX , TSHTRANS I  ENT  SOLUTION  FOR  INSULATED  FLAT  PLATE,  RAH 
IP  FUNCTION  ON  BOUNDING  FLUID! 

30  FORMAT' INI ,TI0,2HK« ,F8.2, !0X,2HC»,F7.4, I 0X,8KDENSlTY»,F8.3l 
AO  FORMAT ( T 1 1 , 1 7HF I LH  COEFFICIENT* ,FI3. 2) 

SO  FORMAT! T 1 1 .6HALPHA* ,F8.A,SX,6HD£L$*  ,F»,3) 

60  FORMAT13FI0. 0,3110) 

TO  FORMAT! 8F 10.0) 

80  FORMAT! 1  HO, T 12, 1 2HE 1  GEN  VALUES) 

30  FORMAT ( T 1 0 , 2HM ! , 1 3 , 3H ) •  ,FtO.S) 

100  FORMAT! IH0.AHTIHE.FI2. 3, SH  SECONDS, I0X, I2HFLUID  TEMP.«.FI2.3> 

110  FORMAT «T2,SF 1 1 .3,SX,5F1 1 .3) 

ISO  FORMAT! El 0 . 1 ! 

READ! 105,10)  ZK.ZC.ZRHO.TF 
WRITE! 106.20) 

WRITE! 108,30)  ZK.ZC.ZRHO 

READ! I0S, 120)  HF1LH 

WRITE! 108. AO)  HFILH 

READ! I  OS. 60)  X.DELX. DELS. KOOE. NODES. N 

READ ( I  OS , 1 0 )  TAU , DELTAU , TAU I . TAU2 

BlOT*HFILH*DELS/( I2.*ZK) 

ALPHA«ZK/!ZC*ZRHO) 

WRITE! 108.50)  ALPHA, DELS 
WRITE! 108.80) 

IFIKODE-I)  2,1.2 

1  READ! I  OS, 70 1  ( AKt I ) , I  *  I ,N) 

GO  TO  3 

2  CALL  E I  GEN (N, AM, BIOT ) 

3  WRITE! 108,90)  ( I , AM! 1 1 . 1  *  I ,N) 

X!«X 

A*TF/TAUI 

6*100. ‘DELS*  *2/ ALPHA 

5  X*XI 

PROD* ALPHA*TAU/3600 . / ! DELS* *2/ 1 AA. ) 

DO  7  J* I .NODES 
SUM*0 , 

DO  6  t*I,N 
ASINI *0SIN!AH! I ! ) 

AS1N2>DSIN!2.*AM( I ) ) 

POWER* AH! I ) • *2*PROD 

IF! POWER. GT. 17A. )  POWER»l7A. 

6  SUM* SUM* ! ASINI • (DEXP! -POWER) *  I )/! AM! 1 )**2*AM( I) *2*ASIN2) I *DCOS!AM( 
II)*X/DELS) 

T! A ) *A* l TAU*B*SUM) 

IF(Tl-J)  *0.01 )  8,8,7 
IF! N- I >  12.12,1 
12  RETURN 
STOP 
END 


.<§ 

! 
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Appendix  B 

LITERARY  SURVEY  OF  THE  NUMERICAL  SOLUTION  OF  THE 
ONE-DIMENSIONAL  HEAT  CONDUCTION  EQUATION 

prepared  by 

University  of  Nevada 
Contract  N60530-67-C-0051 
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PREFACE 


The  bibliography  for  this  suaur y  comprises  a  relatively  small 
part  of  the  complete  bibliography  assembled  under  the  contract.  Only 
the  articles  and  papers  which  appeared  to  deal  directly  with  the 
questions  at  hand,  and  which  were  available,  are  discussed  here.  Many 
of  the  works  were  not  available  at  this  University.  The  remainder  of 
the  bibliography  is  intended  to  offer  a  somewhat  wider  range  of  ref* 
erences  for  information  having  possible  application  in  the  numerical 
treatment  of  heat  conduction  problems. 

The  bibliography  is  arranged  alphabetically  with  a  brief  abstract 
for  each  title.  The  bibliography  is  then  categorized  under  broad 
headings.  Some  titles  may  appear  unuer  more  than  one  heading.  Even 
though  the  contract  requested  a  literature  search  covering  the  area 
of  the  implicit  numerical  solution  of  the  one-dimensional  transient  heat 
transfer  problem  some  explicit  methods  were  considered..  Many  works  com* 
pared  the  methods  so  some  titles  covering  explicit  methods  are  included. 
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NOMENCLATURE 


a  «  space  ind*x  (node  number) 
n  ■  tin*  ind*x  (n  •  1,2,3, 
t  *  time 

at  •  time  increment 

tn  »  elapsed  tine  •  nit  (unifora  at) 

T  ■  temperature 

ax  •  distance  increment 

a  »  thermal  diffusivity  •  L. 

cp 

k  ■  thermal  conductivity 
c  •  specific  heat 
p  *  density 
S  •  Fourier  modulus  • 

L 

N„,  »  Biot  nusber  •  iii. 

k 

L  »  Length 


I 

i 

? 


n«l 

»  t  at.  (arbitrary  at)tQ  »  o 
k»o 


j 

i 

l 
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j  RESULTS  OF  LITERATURE  SEARCH 

i 

PAST  I 

Errors  ia  Inplicit  Finite  Difference  Solutions  of  the  One  Dimensional 
Best  Conduction  Equation 
Definitions 

The  following  definitions  will  be  used,  consistent  with  Anderson  and 
Botje1, 

Roundoff  Error 

This  error  is  caused  by  the  fact  that  all  numbers  used  in  computation 
must  be  rounded  to  a  Manageable  nuaber  of  digits.  The  error  can  become 
significant  after  long  computations  in  which  each  calculation  is  dependent 
on  the  results  of  the  prsviou.  calculation. 

Convergence  Error 

This  is  the  error  caused  by  not  completely  satisfying  the  simultaneous 
equations  when  using  an  iterative  solution.  Increasing  the  number  of 
iterations  decreases  this  error. 

Tine  Truncation  Error 

TMs  arror  can  be  visualised  as  arising  from  the  assumption  that 

temperature  is  a  linear  function  of  time  over  each  time  step;  or  it  can 

be  seen  to  result  from  the  fact  that  certain  high  order  terms  are  neglected 

j  in  the  Taylor  series  type  of  development  of  the  finite  difference  approxi- 

s  mation  to  the  tine  derivative  in  JT  ■a32T.  Reducing  the  size  of  the  time 

i  TT  lx2 

step  reduces  this  error. 

i  Space  Truncation  Error 


Thia  error  can  be  thought  of  as  arising  from  the  assumption  that  temp- 
arature  is  a  linear  function  of  distance  over  each  space  increment;  or  it 
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i 


can  ba  visualized  u  a  consequence  of  the  fact  that  sow  high  order  tens 
are  ne|lected  in  a  Taylor  series  development  of  the  finite  difference 

approxiwtion  to  the  space  derivative  in  -Ji  .  j 

An  additional  type  of  error  is  described  by  Schneider  and  Fox9: 

Error  Caused  by  Boundary,  Discontinuity  ■ 

I 

The  physical  interpretation  of  this  error  is  the  saw  as  that  for 
truncation  error,  given  above,  except  that  near  a  boundary  undergoing  a  1 

step  change  of  temperature,  the  assumption  of  linear  teaperature  variation 
is  poor,  and  the  error  increases  sharply.  It  can  ba  decreased  by  using  a 
suitable  average  between  the  upper  and  lower  step  temperatures  and  by 
reducing  the  length  of  the  tiw  interval.  This  error  decreases  as  tiw 
increases,  in  stable  finite  difference  representations.  In  the  Taylor 
series  development,  this  error  is  seen  to  result  from  the  fact  that  a 

step  temperature  increase  represents  a  discontinuity  in  il*  whereas,  | 

»t  Jjf 

the  Taylor  Series  wthod  assumes  the  derivatives  to  ba  continuous. 

The  foregoing  definitions  apply  to  both  the  explicit  and  implicit 
forms  of  finite  difference  representation. 

SUH3ARY 

Roundoff 

Roundoff  error  is  not  inherent  in  a  finite  difference  approximation, 
but  is  a  type  of  error  associated  with  most  numerical  computations.  It 
can  become  significant  in  calculations  involving  repeated  use  of  rounded 
nuabers.  A  fact  to  be  considered  when  computing  with  finite  differences  is 
that  roundoff  error  is  nearly  always  present,  so  that  even  if  truncation 
error  is  completely  eliminated,  the  solution  cannot  be  entirely  free  from 
error.  The  techniques  to  be  described  for  reducing  truncation  error, 
consequently,  have  no  value  if  the  roundoff  error  is  large  enough  to  over* 
shadow  the  effects  of  such  techniques. 
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Truncation 


Truncation  error  rtpresents  the  difference  between  the  analytical 
solution  and  the  finite  difference  solution  of  a  problem.  Frequently,  a 
Taylor  series  expansion  is  used  to  develop  finite  difference  representations 
and  numerical  estisates  of  truncation  error, 

The  backward  difference  expression  for  developed  from  a  Taylor 
series  expansion,  is  *  *  ’  * 

Hie  finite  difference  expression  for  JT,  derived  by  use  of  a  Taylor 

W 

series  expansion,  is 

f  . 

In  forming  the  simple  backward  difference  implicit  approximation  to  the 
heat  equation,  the  terms  above  which  contain  partial  derivatives  are  neg¬ 
lected.  If  the  neglected  terms  are  considered  to  oe  the  truncation  error, 
the  time  truncation  error  in  the  approximation  to  the  heat  equation  is : 
-At/32T\  ♦  (At)ft3T\  . . 


-  At  /3ZT\  ♦  CAtJWTV  ... 

3atSFJK« 

and  the.  space  truncation  error  is 

(Ax)2/3‘,tV  ♦'•••  E 

lz  \nr)m,n  3<>0  \3x*/m,n 


Evaluation  of  thesa  terms  requires  a 


knowledge  of  the  analytical  solution,  but  if  they  were  to  be  computed  for 
several  cases  having  known  analytical  solutions  and  found  to  have  similar 
values  in  each  case,  their  use  might  be  extended  to  cases  having  no  known 
analytical  solutions.  Kardas11  evaluates  error  terms  for  the  case  of  the 
infinite  plate  with  uniform  initial  temperature  subjected  to  equal  step 
changes  of  surface  temperature  at  time  zero.  The  error  terms  are  found 
using  the  known  analytical  solution  for  the  temperature  distribution,  and 
the  results,  given  as  "error  parameter",  are  plotted  against  NBi>  with  0 
as  parameter.  The  curves  presented,  illustrate  total  truncation  error  only. 
An  example  is  worked  to  illustrate  the  use  of  the  curves,  but  no  indication 
of  the  accuracy  with  which  the  derived  correction  approximates  the  true 
difference  between  analytical  and  numerical  solutions  is  given. 
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Treed  and  Rallis8  have  expanded  the  method  beyond  consideration  of 
the  higher  order  derivative  tents  alone.  Let  •  the  central  difference 
operator  with  respect  to  a  *  T^j^  -  2Ta(I1  ♦  Tm.1<n  and  let  Vn  ■  the  back¬ 
ward  difference  operator  with  respect  to  n  ■  T,^  -  Tmjn„i.  Then,  replacing 

3;T  ■  by  the  backward  difference  implicit  representation,  including 

ax2'  cT  3x 

high  order  derivative  terms: 

6Vn,n  “  Vn  Tm,n  +  Um,n  where  Tm,n  «P™«nts  the  exact  solution 
of  the  differential  equation  and  Uffltn  represents  the  high  order  derivative 
terns  of  the  Taylor  series,  which  are  neglected  in  finite  difference  cal¬ 
culations. 

Also,  consider  f2n  h'B^n  ■  1  ?n  WB^n  which  is  the  difference  equation 
actually  solved  in  finite  difference  calculations.  Let  WB>n  be  the  exact 
solution  of  the  difference  equation;  then  the  truncation  error  Vm<n  is 
defined  to  be  the  difference  between  the  exact  solution  of  the  differential 
equation  and  the  exact  solution  of  the  difference  equation,  or  Vm  n  (trun¬ 
cation  error)  «  Tm<n-WB^n.  This  can  be  given  by  the  difference  between 
the  two  previous  equations  as  Vm>n  -  1_  Vn  Vm<n  ♦  Um>n.  From  this  equa¬ 
tion,  error  estimates  can  be  made  at  each  nodal  point.  The  example  used 
by  Treed  and  Rallis  to  illustrate  the  method  is  the  infinite  plate  at 
uniform  temperature,  subjected  to  identical  step  temperature  changes  at 
the  surfaces.  The  analytical  solution  is  compared  to  the  backward  difference 
implicit  solution,  and  the  errors  predicted  by  the  authors’  method  are 
compared  to  the  errors  predicted  by  computation  of  U  alone.  The  authors’ 
method  yields  better  error  estimates  than  doos  consideration  of  Um  „  only,- 
except  at  the  node  nearest  the  surface.  The  lack  of  improvement  at  that 

point  is  attributed  to  boundary  error.  As  in  the  previous  paper,  time  and 
space  truncation  errors  "re  not  separated,  and  only  total  truncation  error 

is  considered. 

A  different  approach  to  truncation  error  estimation  has  been  used  by 
another  group.  In  this  method,  the  analytical  solutions  for  certain  cases 
are  compared  directly  with  numerical  solutions,  in  an  effort  to  ascertain 


I 

1 

\ 

i 


$ 

% 


5 

s 

3 


I 


4 


1 


74 


_ _ _ _ WC.T?  5143 

truncation  error  relationships  which  can  be  extended  to  cases  for  which 
analytical  solutions  are  not  available.  The  tine  and  space  truncation 
errors  are  considered  separately.  Anderson  and  Botje1  evaluate  space 
truncation  error  for  the  case  of  the  seni  infinite  slab  with  step  temperature 
change  at  the  surface.  Time  truncation  error  is  reduced  to  an  insignificant 
level  by  employing  very  small  time  steps  in  the  numerical  solution.  The 
two  nodes  closest  to  the  surface  are  considered,  and  spatial  truncation 
error  is  shown  to  be  significantly  lower  at  the  second  node  than  at  the 
node  adjacent  to  the  surface.  Curves  art  presented  illustrating  truncation 
error  plotted  against  Fourier  modulus,  using  the  Biot  number  as  a  parameter. 

Anderson,  Slonneger  and  Craybeal'*  study  total  and  time  truncation 
error  for  the  semi-infinite  solid  and  for  the  infinite  plate  with  step 
temperature  change  at  the  surface.  Spatial  truncation  error  is  datermined 
by  following  the  method  used  earlier  by  Anderson  and  Botje1.  Then  other 
numerical  solutions  are  obtained  using  an  arbitrary  time  increment  in  each 
solution.  Subtracting  cne  previously  determined  space  truncation  errors 
from  the  total  truncation  errors  of  these  solutions  yields  the  time  truncation 
errors. 

The  conclusions  of  the  study  arc,  in  summary: 

1.  by  judicious  use  of  the  curves  presented,  values  of  space  and  time 
increments  can  be  chosen  in  such  a  way  that  truncation  error  is 
minimized. 

2,  Truncation  error  curves  for  the  semi-infinite  bodv  and  the  infinite 

plate  closely  approximate  each  other,  indicating  that  the  results 
of  the  paper  can  be  applied  directly  to  other  configurations  with¬ 
out  causing  large  errors. 

3.  Time  truncation  errors  of  significant  magnitude  occur  at  the 
second,  third,  and  fourth  nodes  from  the  surface,  which  contrasts 
with  the  limited  influence  of  space  truncation  errors,  mentioned 
above. 

4,  The  curves  illustrate  the  fact  that  decreasing  the  time  increment 
will  decrease  the  total  truncation  error. 
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This  p*p»r  contains  such  data,  in  the  for*  of  curves,  relating  spatial, 
tie*,  and  total  truncation  error  to  tiaa  and  space  increment,  to  position 
of  a  node,  and  to  elapsed  tiae. 

Space  trunction  error  alone  is  treated  in  the  paper  of  Murray  and 
Landis1'*,  In  their  analysis,  only  tht  space  derivative  of  the  conduction 
equation  is  replaced  by  a  finite  difference  representation,  Khan  the 
resulting  expression  is  applied  to  a  nodal  network,  there  is  obtained  a 
systeN  of  simultaneous  first  order  ordinary  differential  aquations.  Solution 
of  the  systas  by  an  exact  method  givas  a  temperature  distribution  fraa  from 
tiae  truncation  error.  The  aethod  of  solution  employed  by  the  authors  in¬ 
volves  reduction  of  the  differential  equations  to  algtbraic  equations  by 
•mans  of  tho  Laplace  transfora,  with  the  final  solution  obtained  through 
matrix  analysis.  The  case  treated  is  e  slab  with  equal  teaperatura  changes 
at  the  faces,  and  both  step  and  raap  temperature  changes  are  considered. 

The  exact  solution  of  tha  heat  conduction  equation  is  compared  with  the 
exact  solution  of  tha  systam  ot  ordinary  differential  equations  to  daternine 
the  spatial  truncation  error.  The  results  are  plotted  as  truncation  error 
versus  Fourier  modulus,  with  the  nunber  of  nodes  as  parameter.  The  effect 
of  the  convective  film  coefficient  is  not  considered  in  the  exsmple  pre¬ 
sented. 

The  methods  of  the  papers  mentioned  above  might  be  used  to  provide 
estimates  of  truncation  error  for  a  given  set  of  constants  used  in  a  back¬ 
ward  difference  implicit  solution  of  the  conduction  equation,  or  they  might 
be  used  to  select  a  set  of  constants  which  would  minlmiia  the  truncation 
error.  Another  possible  way  to  decrease  truncation  error  is  to  use  one  of 
the  other  implicit  difference  equations,  possessing  lower  inherent  truncation 
error  than  does  the  simple  backward  difference  equation.  Two  of  these, 
the  Crank-Nicolson  and  the  Crandall  equations,  will  be  discussed  briefly. 

The  Crank-Nicolson  equation  is: 

*Tm-l ,n+l  *  2jTm,n*l  *  Tm*l,n«-1  *  Tm-l,n  *(£  -2)Tra,n*Tm*l,n* 
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This  equation  is  discussed  by  Douglas6,  Fox9,  f.aumer10,  and  Campbell, 

Kaplan  and  Moore4.  Douglas  and  Fox  show  mathematically  that  the  truncation 
error  of  this  equation  is  less  than  that  of  the  backward  difference  equation. 
Douglas  shows  the  Crank-Nieoison  equation  to  be  convergent  for  any  value 
of  0.  He  comments  that  the  Crank-Nicolson  equation  provides  considerably 
increased  accuracy  over  the  backward  difference  equation  with  a  small 
increase  in  computation;  but  he  cautions  that  a  lack  of  smoothness  in  the 
solution  of  the  heat  equation  will  retard  convergence  ot  tne  Crank-Nicolson 
solution  to  a  greater  degree  than  convergence  of  the  backward  difference 
equation  will  be  retarded. 

Cauner  compares  the  backward  difference  and  Crank-Nicolson  equations, 
with  regard  to  accuracy  and  stability,  by  applying  both  methods  to  a  prac¬ 
tical  problem  whose  analytical  solution  is  known.  The  problem  considered 
is  the  infinite  plate  with  equal  step  changes  of  temperature  applied  at 
each  face.  The  results  are  presented  as  curves  of  temperature  vs.  time, 
with  the  reciprocal  of  the  Fourier  modulus  as  parameter.  The  curves  illus¬ 
trate  several  points: 

1.  Although  both  numerical  solutions  converge  for  8  »  4  (largest  value 
of  8  used) ,  convergence  is  not  rapid,  and  neither  of  them  is  an 
accurate  approximation  to  the  analytical  solution  at  early  time  for 
such  a  high  modulus. 

2.  Using  8  ^1/4,  both  numerical  solutions  show  rapid  convergence. 

3.  With  8  •  1/4,  the  Crank-Nicolson  equation  offers  slightly  improved 
accuracy  compared  to  the  backward  difference  equation,  but  the 
drawings  do  not  permit  a  precise  evaluation  to  be  made. 

The  Crandall  equation  is  discussed  by  Crandall9,  Douglas6,  and  Campbell, 
Kaplan  and  Moore*.  It  is  related  to  the  Crank-Nicolson  equation,  differing 
only  in  that  different  constants  are  used.  The  Crandall  equation  can  be 
expressed  as  follows: 

Tm,n*l  *  Tm,n  *  1/?2  "  ^6)  ^Tm*l,n»l  "  2Ta,n*l  *  Tm-l,n*l^ 

♦  1/2  (8  ♦  1/6)  (Vlfn  -  2Tm,n  ♦  Vl>n) 
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Cam pbell,  Kaplan  and  Moor*  compare  th*  Crandall  and  Crank»Hicolson 
equations  with  th*  analytical  solution  for  th*  temperature  distribution  in 
an  infinite  plate  subjected  to  aqual  step  temperature  changes  at  the  surfaces. 
Results  are  presented  for  only  on*  value  of  Fourier  modulus,  which  is  given 
in  th*  body  of  the  paper  as  1/2,  but  is  reported  as  2  in  the  conclusions 
section.  The  Crandall  equation  is  shown  to  offer  increased  accuracy  con* 
pared  to  the  Crank-Nicolson  aquation,  but  no  comparison  of  computation 
time  is  given.  Another  fact  demonstrated  by  Campbell,  Kaplan  and  Moore  is 
that,  as  the  time  and  space  increments  are  decreased  (9  maintained  constant), 
tue  error  of  the  Crandall  equation  is  reduced  more  rapidly  than  that  of 
the  Crank-Nicolson  equation. 

As  mentioned  above,  Campbell,  et.al.,  present  truncation  error  for 
only  one  value  of  Fourier  modulus.  Generally,  truncation  error  is  different 
for  different  values  of  Tourier  modulus,  and  it  is  of  interest  to  note  that 
Crandall2  recommends  the  us*  of  9  »  0.2236  in  the  Crandall  difference 
equation  as  the  value  which  should  give  the  smallest  truncation  error. 

Another  approach  to  reducing  the  truncation  error  in  finite  difference 
approximation  is  an  extension  of  the  process  known  as  Richardson's  deferred 
approach  to  the  limit,  discussed  by  Douglas,  Fox  and  Liebmann.  Douglas6, 
and  Fox9,  demonstrate  the  mathematical  validity  of  the  method,  and  Lieb¬ 
mann14,  illustrates  its  application  to  a  practical  problem.  In  the 
example  given  by  Liebmann,  the  Richardson  technique  is  used  to  reduce 
time  truncation  error  in  the  backward  difference  equation.  The  procedure 
will  be  described  here:  The  backward  difference  solution  is  carried  through 
using  a  given  time  increment,  At.  Then  the  problem  is  solved  again  using 
the  same  equations,  but  with  the  time  increment  doubled  in  length  (2At). 

Then,  for  a  given  space  node  at  a  given  time,  the  temperature  T(At)  obtained 
by  using  At,  is  corrected  by  adding  to  it  the  difference  between  T(Atj 
and  T(2At) ,  where  T(2At)  is  the  temperature  obtained  using  the  time  in- 
crement(2At) .  Thus  T(corrected)  ■  T  (At)  ♦  [T  (At)  -  T  (2At)].  Douglas6 
mentions  that  the  technique  can  be  used  to  reduce  spatial  truncation 
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error  rather  than  time  truncation  error,  by  using  a  linear  coatoination  of 
two  solutions,  each  of  which  uses  a  different  value  of  space  increment  Ax, 
with  At  being  the  same  in  both  solutions.  He  also  comments  that  both 
spatial  and  time  truncation  error  can  be  reduced  by  using  a  linear 
combination  of  three  solutions,  and,  furthermore,  that  the  basic  idea  can 
be  applied  to  equations  other  than  the  backward  difference  equation. 

Boundary  Error 

Boundary  error  is  introduced  at  the  start  of  calculations,  end  it 

decreases  as  the  solution  is  advanced  in  the  time  dimension,  provided  the 

finite  difference  equation  being  used  is  stable.  Implicit  finite  difference 

equations  are  stable  for  all  values  of  the  Fourier  Modulus;  however,  the 

boundary  error  is  often  excessive  at  early  ties  steps  unless  the  Fourier 

modulus  is  snail.  This  requires  that  tha  ratio  of  the  tine  step  to  the 

distance  step  be  relatively  snail,  which  naans  that  a  large  nunber  of 

calculations  is  required  to  cover  a  given  tine  interval.  As  tine  passes, 

the  requirenent  for  snail  At  it  diminished  because  of  the  inherent 

lx 

decrease  in  boundary  error.  To  reduce  the  amount  of  computation,  the  tine 
increnr  *.  can,  therefore,  be  increased  as  the  solution  progresses,  thus 
increasing  0  also. 

Douglas6  and  Douglas  and  Gallia7  discuss  variable  tine  steps,  and 

present  two  schemes  for  increasing  the  length  of  the  tine  step  in  a  systems* 

tic  way.  The  first  method  results  in  a  linear  increase  of  the  tins  incre* 

nent  as  the  solution  of  tho  difference  equation  advances  in  the  tine  din* 

ension.  To  determine  the  tine  interval  to  bo  used  between  tine-tn  and 

tine  tRtl,  use  Atn  »  (a  ♦  btn)  (Ax)2  whore  a>o  and  b>o  and  tn(n»l,2,3,) 
n-1 

■  £  Atk,  o<k<«,  (t  a  o).  The  second  method  causes  the  length  of  the 
k»o 

tine  step  to  grow  exponentially.  The  tin*  interval  for  use  between  tR 
and  tjj^j  is  given  by  Atn  a  (Ax)2  e  *tn  where  o<a<^2  and  tn  is  given  by  the 
* unnation  above.  Douglas  and  Gallie  discuss  the  use  of  these  relations 
as  to  the  backward  diffsrence  implicit  equation,  and  show  that  the 

rate  of  convergence  of  the  solution  is  not  reduced  by  their  use,  which 
iaplies  that  the  accuracy  at  a  given  elapsed  tint  is  not  reduced. 
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PART  II 


Some  Method*  of  Solvinf  Implicit  Finis*  Diff*r*nc*  R*pr*s*nt*tions 
of  th*  Heat  Conduction  Equation. 


Two  methods  employed  for  solving  th*  systems  of  aquations  resulting 

fro*  us*  of  th*  implicit  difference  equations  ar*  iteration  and  Gaussian 
elimination.  Two  papers  comparing  different  variations  of  thwse  processes 

and  so**  references  containing  special  techniques  for  use  in  solving  the 
systew  will  be  Mtiticned. 

2 

Anderson,  Botje,  and  Koffel  enploy  the  backward  difference  equation  in 
a  coaputer  program,  using  Gauss-Seidel  iteration  to  solve  the  simultaneous 
equations,  for  a  two  diMnsional  network  of  as  nany  as  200  nodes.  Two 
sch*Ms  are  used  in  coebination  to  accelerate  convergence  of  th*  iteration, 
and  they  will  be  described  later.  Gaussian  elimination  was  initially 
considered  by  these  authors  for  solving  the  simultaneous  equations,  but  it 
was  found  that  excessive  computation  tin*  would  be  required  for  such  a 
large  r.uaber  of  equations.  The  authors  describe  thoroughly  th*  development 
and  application  of  this  prograa  and  give  results  for  several  industrial 
problems  it  has  solved. 

Several  methods  are  available  for  accelerating  the  convergence  of 

iterative  processes,  thus  reducing  computation  time,  Anderson,  Botje  and 

Koffel2  discuss  their  experience  with  a  combination  of  two  such  devices  in 

connection  with  an  implicit  heat  transfer  program,  commenting  that  time 

savings  of  as  much  as  75  percent  have  been  obtained.  The  first  of  these 

schemes  involves  extrapolation,  the  initial  value  for  an  iteration  process 

being  obtained  by  extrapolating  fron  the  results  of  the  two  proceeding 

iteration  steps.  The  second  method  is  an  adaptation  of  the  technique 

described  by  Wegstein17,  which  is  very  similar  to  the  Aitken  f> 2  process. 12>15 

In  the  Wegstein  method,  the  value  of  the  unknown,  x  at  the  (k*l)  iteration 

step,  is  corrected  as  follows:  X.  .  (corrected)  •  X.  .  -  (xk*l  *  fy)2 

xk*l  *  *h*  xk-l 


80 


The  Kegstein  technique  accelerates  convergence  for  sisple  iterative  processes 
which  exhibit  either  oscillatory  or  monotonic  convergence,  and  alto  can  be 
used  to  force  the  convergence  of  simple  iterative  processes  which  show 
oscillatory  or  monetonic  divergence,  Anderson,  Botje  and  Koffel  stake  two 
consents  about  this  technique  based  on  their  experience  with  its  use  in 
their  heat  transfer  computer  progra at: 

1.  If  this  acceleration  adieus  is  applied  as  often  as  every  third 
iteration  sweep,  the  extra  stachine  tise  required  to  coiyute  the 
correction  stay  exceed  the  tine  saving  accomplished  by  the  acceler¬ 
ation. 

2,  A  satisfactory  swthod  for  determining  the  number  of  iteration 
sweeps  between  applications  of  the  acceleration  correction  is  to 
set  the  number  of  sweeps  between  corrections  equal  to  the  nuaber 
of  nodes  froa  the  boundary  to  the  deepest  node  in  the  sys tea. 

A  third  acceleration  technique  is  attributed  to  Steffensen.  It  is 
described  briefly  by  Prager15,  who  states  that  one  application  of  this 
method  has  the  same  effect  as  three  successive  applications  of  Aitken's  42 
process. 

A  modified  fora  of  Gaussian  elimination,  suitable  for  solving  the 
tridiagonal  system  of  equations  obtained  from  application  of  the  backward 
difference  inplicit  equation,  is  described  by  Douglas6,  and  Wilkes**.  It 
is  based  on  the  fact  that  the  difference  equation  at  each  interior  node 
contains  three  unknowns,  and  those  at  the  boundaries  contain  two  unknowns. 
Using  this  property,  a  general  expression  is  developed  relating  the 
unknown  temperature  at  each  node  to  the  temperatures  of  adjacent  nodes. 

An  expression  is  obtained  giving  the  temperature  of  the  final  node  of  the 
system  explicitly.  Then  the  temperature  of  each  node  in  turn  is  calculated, 
beginning  at  the  final  node  and  proceeding  toward  the  first  node. 
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PART  III 


Truncation  Error  of  Runge-Kutta  Methods 

If  the  second  order  (spstisl)  derivative  alone  is  replaced  by  a  finite 
difference  expression,  the  one  dimensional  heat  conduction  equation  becomes: 

£L  •  «  (T«*l)n  -  2T  ♦  T^j  n).  When  this  equation  is  applied  to  a 

dt  (dx)z  * 

system  of  heat  transfer  nodes,  a  system  of  simultaneous  ordinary  differential 
equations  is  obtained.  Solution  of  the  simultaneous  system  by  an  exact  method 
would  yield  temperature  values  free  from  time  truncation  error.  If  the 
fystem  is  solved  by  a  numerical  scheme,  truncation  error  will  be  introduced, 
because  such  methods  are  based  on  approximate  relationships.  For  example, 
the  Runge-Kutta  method,  a  widely  known  device  for  solving  ordinary  differential 
equations,  is  developed  from  a  Taylor  series  expansion,  and  contains  a  trun¬ 
cation  error  because  of  the  fact  that  high  order  terms  of  the  series  are 
ignored.  An  estimate  of  the  truncation  error  of  one  set  of  Runge-Kutta  formulas 
is  given  by  Prager*5. 

To  solve  ■  f(u,v)  with  the  initial  condition  v(uQ)  ■  vQ,  let  h  ■ 
interval  length,  p  -  number  of  steps,  V  «  approximation  to  v  given  by  the 
Runge-Kutta  equation;  the  most  widely  used  Runge-Kutta  equation  is 

vp»l  ■  vp  ♦  l/6(k1  ♦  2k2  ♦  2k3  ♦  k4)  where 
kj  .  hf(up,Vp) ;  k2  -  hf  (Up  ♦  h,  Vp  ♦  *1) 
k3  -  hf(up  eh,  Vp  ♦  Jj2);  k4  -  hf(U?  ♦  h,  Vp  ♦  k3) 

The  truncation  error  is  estimated  to  be  -1  (Vp+3„Vp)  and  the  corrected  value 
of  vpei  is  vp»l  *  yj  Cvpel  *  V  ♦ 

Lance  12 ,  describes  a  modified  Runge-Kutta  procedure  for  digital  computers 
which  automatically  adjusts  the  interval  length  to  maintain  a  predetermined 
truncation  error. 

Vl  •  Vp  ♦  7<kl  ♦  4k4  ♦  k5) 

k,> 

kj  •  Jhf  (Up  ♦  Vp  ♦  h  ♦  ^2 


k,  -  i.  hf  (U 


3  r  T- 
♦  v„  ♦  ♦ 


p  -g- 


k  -  i  hf  (u  ♦  h,  V_  ♦  3kl  -  9*3  ♦  6k,) 
5  3  p  p  -ip-  4' 
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The  truncation  error  c  is  estimated  by 

Sc  ■  kj  »  9k3  ♦  4k.  ■  k5,  If  tha  right  hand  sida  of  this  aquation 

T*  4  T 

is  graatar  than  fiva  times  tha  allowable  err or,  h  is  halved  and  tha  computation 
for  tha  step  is  repeated;  but  if  tha  right  sida  is  lass  than  S/32  of  the 
allowable  error,  tha  interval  nay  be  doubled  and  the  coaputation  repeated. 
Although  this  process  would  seen  to  require  considerable  extra  tine,  Lance 
reports  that  it  usually  reduces  the  co^iutation  tine  required  to  attain  a 
given  accuracy  of  solution  by  about  20  per  cant.  He  accounts  for  the  re¬ 
duction  by  mentioning  that  when  a  fixed  interval  length  is  used  throughout 
the  solution,  this  length  is  usually  deliberately  underastiaated  to  insure 
accuracy,  thus  requiring  the  use  of  more  steps  than  necessary.  The  self 
adjusting  procedure  eliminates  more  than  enough  of  these  extra  steps  to 
offset  the  time  required  for  the  extra  numerical  manipulations  it  requires, 
thus  effecting  an  overall  saving  of  computation  time. 
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Vol,  13,  part  3,  Aug.  1960,  p.  374-84 

This  article  treats  the  use  of  automatic  computers  for  the  numerical 
solution  of  the  cylindrical  heat  conduction  problem.  It  is  shown 
that  accurate  solutions  can  be  obtained  easily  and  rapidly  using  the 
Crank-Nicolson  implicit  method.  Attention  is  given  to  the  adequacy  of 
the  finite-difference  representation  near  a  singularity  at  the  boundary. 
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i;  the  Heat  Conduction  Equation,"  The  Quarterly  Journal  of  Mechanics 
and  Applied  Mathematics,  Vol,  4,  p,  209-22,  19S1 

The  equation  considered  is  3u/3t  •  k32v/3x2.  The  authors  make  the 
transformation  v  ■  3u/3t  ♦  k  a2u/nx2 ,  which  gives  the  equation 
32u/3t2  -  k23'*u  /3x“*  »  0,  The  boundary  conditions  are  also  transformed 
and  new  ones  added. 

5.  Allione,  M.0, ,  "Comparative  Study  of  Runge-Kutta  and  Lanczos  Numerical 
Integration  Methods",  Rept.  No.  U2421,  Contract  AF  19  628  562,  ESD 
TDR63  662,  Aeronutronic,  Newport  Beach,  Calif.,  9  Jan.  1964,  24  p, 

AD-429  958 

The  two  methods  are  compared  in  solving  the  system  of  ordinary 
differential  equations  associated  with  the  Variation  of  Parasenters 
formulation.  Results  of  ephemeris  calculations  using  each  method 
are  compared  with  a  standard  to  determine  the  relative  error  growth. 
Conclusions  are  drawn  regarding  the  relative  merits  of  the  two  methods. 

6.  Anderson,  J.T.,  "Review  of  Digital  Computer  Heat  Transfer  Programs" 

ASME  Paper  6S-WA/HT-48,  7  p. 

This  paper  is  a  review  of  the  available  steady  state  and  transient 
programs.  Emphasis  is  given  to  the  capabilities  and  limitations  of 
general  purpose  programs,  both  explicit  and  implicit.  Indications  of 
computer  time  are  given  and  discussion  of  the  error  magnitude  is 
included. 

7.  Anderson,  J.  T.,  and  Botje,  J.M.,  "Spatial  Truncation  Error  Analysis, 
"ASME  paper  No.  62-HT-27 

A  method  is  presented  for  evaluating  spatial  truncation  errors  in  a 
finite  difference  solution  of  a  parabolic  partial  differential  equation. 
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A  one -dinars*  A  on  transient  hsst  transfer  problem  is  used  m  an  example. 
Curves  are  presented  for  rapid  evaluation  of  the  spatial  truncation  error. 

8.  Anderson,  J.T.,  Botje,  J.M.,  Koffel,  W.K.,  "Digital  Computer  Solution 
of  Co^lex  Transient  Heat  Transfer  Problems,"  W.Va.  Univ.  Bulletin, 

Engr,  Experiment  Station,  Technical  Bulletin  No.  62,  26  p. 

The  authors  describe  a  comprehensive  computer  program  for  heat  transfer 
transient  problems  involving  convection,  conduction,  contact  resistance, 
solid  and  gaseous  radiation,  surface  flux  and  internal  heat  gsneration. 

The  Liebman  backward  time  step  approximation  was  used  in  developing 
the  program,  and  the  difference  equations  obtained  by  the  Liebman  method 
are  presented  with  a  discussion  of  methods  used  to  solve  them  on  a 
digital  computer. 

9.  Anderson,  J.T.,  Slonnegar,  R.D.,  Creybeal,  C..E,,  'Truncation  Error 
Analysis  for  Transient  Heat  Transfer  Calculations."  unpublished 

"A  total  and  time  truncation  error  analysis  of  numeric  solutions  of 
parabolic  partial  differential  equations  is  herein  reported.  Math¬ 
ematical  models  include  transient  heat  conduction  in  an  infinite  body 
and  an  infinite  plate  with  (1)  convective  heat  transfer  and  (2)  a  step 
function  temperature  change  on  the  surfaces.  Results  are  tabulated  as 
well  as  shown  graphically."  author's  abstract 

10.  Bsrakat,  H.Z.,  and  Clark,  J.A.,  "On  the  Solution  of  the  Diffusion 
Equations  by  Numerical  Methods,"  ASME  Journal  of  Heat  Transfer,  Vol.88, 
p.  421-27,  1966 

Author's  introduction:  "An  explicit  finite  difference  approximation 
procedure  which  is  unconditionally  stable  for  the  solution  of  the 
general  multi-dimensional,  non-homogeneous  diffusion  equation  is 
presented.  This  method  possesses  the  advantages  of  the  implicit  methods, 
i.e.,  no  severe  limitation  on  the  size  of  time  increment.  Also  it  has 
the  simplicity  of  the  explicit  methods  and  employs  the  same  "marching" 
type  technique  of  solution.  Results  obtained  by  this  method  are 
compared  with  the  exact  solution  and  with  those  obtained  by  other  finite 
difference  methods.  For  the  examples  solved  the  numerical  results 
obtained  by  the  present  method  are  in  closer  agreement  with  the  exact 
solution  than  are  those  obtained  by  other  methods" 

11.  Bellman,  R.  Kalaba,  R.,  Xotkin,  B.,  "On  a  New  approach  to  the 
Computational  Solution  of  Partial  Differential  Equations,"  Proceedings 
of  the  National  Academy  of  Sciences  of  the  USA,  Vol.  48,  P.  1325-27,  1962 

This  article  discusses  a  modification  of  the  usual  finite  difference 
approach  to  the  numerical  solution  of  partial  differential  equations. 

The  idea  is  that  the  computational  algorithm  should  exhibit  as  closely 
as  possible  the  properties  of  the  actual  solution;  for  example,  if  the 
actual  solution  is  bounded  and  non-negative,  this  should  be  evident  from 
the  algorithm.  The  method  is  illustrated  by  a  problem,  and  the  numerical 
results  are  discussed. 

12.  Buglia,  J.J.,  and  Brinkworth,  H.,  "A  Comparison  of  Two  Methods  for 
Calculating  Transient  Temperatures  for  Thick  Walls,"  NACA  Technical 
Note  4343,  ID  p.,  Aug.  1958 

This  paper  compares  Hill's  method  (NACA  Tech.  Note  4105)  with  Dusinberre's 
method.  In  Hill's  method,  finite  differences  are  taken  only  in  the 
time  variable,  the  equations  used  being  already  integrated  with  respect 
to  distance.  The  authors  conclude  that  Hill's  method  is,  practically, 
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an  exact  method  and  is  faster  than  Dusinberre's,  only  the  two  surface 
temperatures  are  needed.  If  the  temperature  distribution  is  needed. 

Hill's  method  is  slower,  but  very  accurate. 

13.  Butler,  R,  Kerr,  E.,  "An  Introduction  to  Numerical  Methods,"  Pitman 
Publishing  Ccrp.,  New  York,  386  p.,  1962 

This  is  an  elementary  text  covering  the  solution  of  algebraic  equations, 
finite  differences,  interpolation,  numerical  differentiation  and  Integra* 
tion,  and  the  solution  of  ordinary  differential  equations. 

14.  Campbell,  B.C,  "An  Investigation  of  the  Accuracy  of  Numerical  Solution's 
in  the  Diffusion  Equation  for  Transient  Heat  Transfer,"  Master's 
thesis,  AFIT  CSF/Phys/64-1,  A.F,  Inst,  of  Technology,  Wright-Patterson 
AFB,  Ohio,  Aug.  1964,  118  p.  AD-60S-489 

This  paper  gives  the  results  of  a  comparison  of  the  accuracy  obtained 
using  two  implicit  finite  difference  representations  of  the  transient 
heat  conduction  equation  in  obtaining  solutions  for  two  basic  heat 
conduction  problems.  The  representations  compared  are  the  Crank-Nicolson 
form  and  a  theoretically  more  optimum  form  called  the  Optimum  Implicit 
form.  It  is  demonstrated  that  the  major  source  of  error  occurs  at  the 
initial  time/space  node  of  the  difference  mesh.  The  Optimum  Implicit 
form  decreases  this  error  and  is  shown  to  be  as  accurate  as  the  Crank- 
Nicolson  form  in  one  problem  and  considerably  more  accurate  in  the  other. 

15.  Carr,  J.W.,  III,  "Error  Bounds  for  the  Rungc-Kutta  Single  Step  Integration 
Process,"  Journal  of  the  Association  for  Computing  Machinery ,Vol.  5, 

p.  39-44,  1958 

This  article  presents  a  mathematical  theorem  for  determining  the  error 
at  each  step  in  a  fourth  order  Runge-Kutta  procedure. 

16.  Collatt,  L.,  "The  Numerical  Treatment  of  Differential  Equations,"  third 
edition,  Springer-Verlag,  Berlin,  554  p.,  1960 

This  is  a  rather  comprehensive  book  with  over  200  pages  devoted  to 
partial  differential  equations,  including  discussions  of  error  propagation, 
node  spacing  and  iterative  methods.  Emphasis  is  placed  on  manual  solution 
of  problems  rather  than  computer  solution. 

17.  Compton,  W.R.,  "An  Extension  of  Present  Numerical  Solutions  for 
Transient  Heat  Conduction,"  NOTS  TP 3361  NAVWEPS  8419,  NOTS  China  Lake, 
Calif.,  Feb.  1964,  26  p.  AD-431  791 

A  fourth  order  technique  for  the  numerical  solution  of  transient  heat 
transfer  equations  involving  conduction,  convection,  and  radiation  is 
presented.  Approximating  parabolas  and  Taylor  series  expansions  are 
used  to  facilitate  the  use  of  fourth  order  difference  equations  and 
Runge-Kutta  techniques. 

18.  Crandall,  S.H.,  "On  a  Stability  Criterion  for  Partial  Difference 
Equations,"  Journal  of  Mathematics  and  Physics,  V,  32,  p.  80-81, 

1953 

The  author  exhibits  a  partial  difference  equation  which  is  unstable  but 
has  stable  characteristics  locally.  He  points  out  that  a  stability 
analysis  based  on  the  procedure  of  O'Brien,  Hyman,  and  Kaplan  (J.  Math, 
and  Phys.  V.  29,  p,  223-51,  1951)  is  not  valid  for  this  type  of  equation. 
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19.  Crank*  J.,  and  Nicolson,  P«,  "A  Practical  Method  for  Numerical  Evaluation 
of  Solutions  of  Partial  Differential  Equations  of  the  Heat  Conduction 
Type,"  Proceedings  of  the  Caabridge  Philosophical  Society,  Vol.  43* 

p.  50*67 *  1247 

The  article  presents  a  comparison  of  three  methods  for  the  solution 
of  the  non-linear  equation  of  heat  flow  in  a  medium  where  heat  is  being 
generated.  The  first  method  consists  in  a  reduction  to  a  system  of 
ordinary  differential  equations  by  approximating  the  time  derivatives 
with  differences.  The  second  method  is  the  same  except  that  the  space 
derivatives  are  approximated  by  differences  rather  than  the  time  deri« 
vatives.  In  the  third  method,  all  derivatives  are  approximated  by 
differences,  and  the  authors  conclude  this  to  be  much  faster  and  more 
satisfactory.  The  third  method  gives  a  system  of  non-linear  algebraic 
equations,  which  is  solved  by  a  combination  iterative  and  step  by  step 
method,  and  several  variations  of  this  method  are  discussed. 

20.  Cudahy,  G.  F„,  "Investigation  of  Accelerating  the  Convergence  of  an 
Implicit  Numerical  Solution  of  Transient  Heat  Transfer  Problems," 

Master's  thesis,  Rept.  No.  GA/PH/65-4-A,  A.F.  Institute  of  Technology, 
Wright-Patterson  AFB,  Ohio,  Aug.  1965,  82  p,  AD-621  273. 

This  paper  presents  the  results  of  an  investigation  of  two  methods  for 
increasing  the  convergence  rate  for  the  two  dimensional,  five  point 
implicit  finite  difference  representation  of  the  diffusion  equation 
of  transient  heat  transfer,  the  methods  being  the  adapted  Kegstein 
technique  and  successive  overrelaxation.  Various  mesh  scanning  techniques 
are  investigated.  An  example  problem  is  used  to  show  that  successive 
overrelaxation  with  a  technique  of  repeatedly  scanning  all  boundary 
values  into  the  finite  difference  mesh  is  the  fastest  method  of  solution 
for  the  equations  used  in  this  study.  Solution  of  28  other  problems 
by  this  method  shows  an  approximate  saving  of  i/3  in  the  number  of 
iterations  over  successive  overrelaxation  with  a  conventional  repetitive 
scanning  technique.  An  "a  priori"  relaxation  factor  related  to  the 
maximum  temperature  gradients  is  obtained, 

21.  Descloux,  J,,  "Note  on  the  Round-off  Errors  in  Iterative  Processes," 
Mathematics  of  Computation,  Vol.  17,  p,  18-37,  1963 

Let  Xj,,i  «  X,j  ♦  F  (Xn)  be  a  scalar  iterative  converging  process.  When 
Xn  is  close  to  the  limit,  F  (Xn)  is  small  and  can  perhaps  be  obtained 
with  higher  absolute  precision  than  X^  Then  the  addition  Xn  ♦  F  (Xn) 
involves  a  round  off  operation.  The  author  shows  that,  for  a  fixed-point 
computer,  an  appropriate  rounding  method  can  impove  the  accuracy  of 
solution.  Appendix  1  gives  analogous  results  for  a  floating-point 
computer.  Appendix  II  deals  with  Aitken's  62  process, 

22.  Douglas,  Jr.,  J,  "A  Note  on  the  Alternating  Direction  Implicit  Method  for 
the  Numerical  Solution  of  Heat  Flow  Problems,"  Proceedings  of  the  American 
Mathematical  Society,  Vol,  8,  p.  409-12,  19S7 

As  originally  presented,  the  alternating  direction  method  of  Peaceman 
and  Rachford  applies  to  rectangular  regions.  The  author  extends  it  to 
regions  with  polygonal  boundaries  where  each  segment  of  the  boundary  is 
parallel  to  one  of  the  coordinate  axes. 

23.  Douglas,  Jr.,  J.,  "A  Survey  of  Numerical  Methods  for  Parabolic  Differential 
Equations,"  Advances  in  computers,  Vol.  2,  p  1-54,  1961,  Academic  Press, 
New  York. 
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Quoted  from  the  euthor's  introduction:  "The  purpose  of  this  survey  is  to 
introduce  a  theoretically  minded,  but  not  highly  mathematically  trained, 
scientist  to  finite  difference  Mthods  for  approximating  the  solutions 
of  partial  differential  equations  of  parabolic  type,"  A  partial  list  of 
topics  is:  Explicit  Difference  Equations,  The  Backward  Difference 
Equation,  The  Crank-Nicolson  Equations,  Comparison  of  Calculation 
Requirements,  Alternating  Direction  Methods. .Higher  Order  Correct  Difference 
Equations 

24.  Douglas,  Jr.,  J.,  and  Gallic,  Jr.,  T.M.,  "Variable  Time  Steps  in  the 
Solution  of  the  Heat  Flow  Equation  by  a  Difference  Equation,"  Proceedings 
of  the  American  Mathematical  Society,  Vol.  6,  p.  787-93,  1955, 

The  authors  consider  the  numerical  solution  of  Ut  *  Uxx,  using  a  back¬ 
ward  difference  equation  known  to  be  stable  for  all  Ax  and  At.  Two  cases 
of  variable  At  are  considered. 

25.  Douglas,  Jr.,  J.,  "The  Solution  of  the  Diffusion  Equation  by  a  High 
Order  Correct  Difference  Equation,"  Journal  of  Mathematics  and  Physics, 
Vol.  35,  p.  145-51,  1956 

The  author  proposes  a  six-point  implicit  difference  scheme  for  the 
solution  of  Ut  ■  Uxx  with  a  smaller  truncation  error  than  the  Crank- 
Nicolson  form.  Because  the  error  is  smaller,  the  new  method  allows 
larger  AX  and  At  to  be  used,  (The  method  is  described  also  in  the  book 
"Advances  in  Computers,  Vol.  II,  p.  26.) 

26.  Douglas,  Jr.,  J.,  "The  Effect  of  Round-Off  Error  in  the  Numerical 
Solution  of  the  Heat  Equation,"  Association  for  Computing  Machinery 
Journal,  Vol.  6,  No,  1,  Jan,  1955,  p,  48-58 

The  artical  presents  an  analysis  of  the  approximation  to  the  heat 
equation  by  the  backward  difference  equation  when  boundary  value  problems 
are  approximated  by  finite  difference  problems.  The  dependence  on  the 
method  of  solving  the  tri-diagonal  equations  is  shown.  It  is  shown  that 
if  linear  equations  are  solved  by  a  normalized  form  of  Gaussian  elim¬ 
ination,  the  procedure  is  stable  against  round-off  error. 

27.  Dusinberre,  G.M. ,  "A  Note  on  the  'Implicit*  Method  for  Finite-Difference 
Heat-Transfer  Calculations,"  ASME  Journal  of  Heat  Transfer,  Vol.  83,  p.  94 
1961 

"The  apparent  advantage  of  the  implicit  method  lies  in  the  possibility 
of  using  relatively  large  time  intervals.  But  this  may  be  accompanied  by 
(1)  considerable  sacrifice  in  accuracy  and  (2)  no  corresponding  saving 
in  digital  time."  This  is  the  author's  introduction  to  the  article,  which 
discusses  points  (1)  and  (2)  above. 

28.  Dusinberre,  G.M.,  "Numerical  Methods  for  Transient  Heat  Flow,"  ASME 
Transactions,  Vol.  67,  p.  703-10,  1945 

A  modulus  is  developed  by  choice  of  which  a  solution  may  be  reached  most 
rapidly  or  alternatively  reached  more  slowly  but  with  greater  precision. 
Criteria  are  developed  for  choosing  the  modulus  to  insure  convergence. 

A  method  is  developed  for  handling  thermal  conductivity  and  heat  capacity 
when  they  vary  independently  with  temperature. 

29.  Elliot,  D.,  "A  Method  for  the  Numerical  Integration  of  the  One-Dimensional 
Heat  Equation  Using  Chebyshev  Series,"  Proceedings  of  the  Cambridge 
Philosophical  Society,  Vol.  57,  p.  823-32,  1961 
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The  aquation  used  is  8S/3t  - 32e/»x 2  (-l<X<l;t>0  with  general  linear 
boundary  conditions  along  x  *  <1.  38/3t  is  replaced  by  a  finite  difference 

approximation  and  the  resulting  system  of  ordinary  differential  equations 
is  solved  by  Clans haw's  method  of  0  in  teres  of  Chebyshev  polynoaials. 

Two  examples  are  worked  and  the  results  compared  with  exact  results. 

The  author  concludes  that  the  present  method  requires  less  computation 
than  the  usual  finite  difference  methods,  but  is  less  versatile  and  not 
so  well  suited  for  coaplicated  equations . 

30.  Elrod,  Jr.,  H.G.,  "New  Finite-Difference  Technique  for  Solution  of  the 
Heat-Conduction  Equation,  Especially  Near  Surfaces  with  Convective  Heat 
Transfer,"  ASME  Trensactions,  Vol.  79,  p.  1S19-25,  19S7 

The  success  of  most  finite  difference  methods  for  transient  heat  conduction 
depends  on  the  existence  of  a  certain  degree  of  uniformity  of  behavior 
of  the  temperature  over  the  time  and  space  intervals  selected  for  coaputation. 
This  often  requires  the  use  of  inconveniently  short  time  intervals. 

This  paper  represents  an  effort  to  develop  a  finite  difference  method 
not  possessing  such  a  defect. 

31.  Eaacns,  H.  to',,  "The  Numerical  Solution  of  Heat  Conduction  Problems, " 

ASME  Transactions,  Vol.  65,  p.  607-12,  1943 

The  author  discusses  the  application  of  the  Southwell  relaxation  method 
to  two  and  three  dimension  steady  state  heat  conduction.  A  transient 
problem  is  also  briefly  considered.  A  short  discussion  on  "Derivation 
of  Difference  Equations  from  Differential  Equations,"  is  included. 

32.  Emmons,  H.  W.,  "The  Numerical  Solution  of  Partial  Differential  Equations" 
Quarterly  of  Applied  Mathematics,  Vol.  2,  p.  173-95,  1944 

The  author  presents  a  detailed  expository  treatment  of  Southwell's 
relaxation  process.  Examples  illmstrste  the  application  to  the  solution 
of  boundary  value  problems  for  the  Laplace,  Poisson  and  other  equations, 

33.  Enig,  J.W, t  "A  Method  for  the  Rapid  Numerical  Solution  of  the 
Heat  Conduction  Equation  for  Composite  Slabs,"  NAVORD  Rept.  6666 
Naval  Ordnance  Lab.,  White  Oak,  Md.,  20  Aug  59,  22  p.  PB  144  193 

Two  boundary  conditions  encountered  in  heating  one  dimensional  double 
slabs  are:  heat  flux  given  at  inner  surface  and  (a)  heat  flux  or 
(b)  temperature  given  at  the  outer  surface.  A  method  is  developed  which 
permits  rapid  calculation  of  (1)  any  interior  point  temperature  and 
(2)  the  outer  surface  temperature  or  flux  for  (a)  or  (b)  respectively, 
without  computing  other  interior  point  temperatures.  The  partial  diff¬ 
erential  equations  are  integrated  in  terms  of  an  arbitrary  outer  surface 
temperature  or  flux  by  a  simple  numerical  scheme.  The  method  performs 
an  exact  integration  over  the  space  dimension,  so  once  the  outer  surface 
temperature  is  determined  the  interior  temperatures  are  computed  by  exact 
formulae.  The  numerical  solutions  are  compared  to  the  exact  solution 
for  accuracy  and  to  other  numerical  schemes  for  speed. 

34.  Forsythe,  G.E.,  and  Wasow,  W.R.,  "Finite-Difference  Methods  for  Partial 
Differential  Equations",  John  Wiley  end  Sons,  Inc.,  New  York-London 
444  p.,  1960 

The  heat  equation  Ut  ■  U  is  considered  (in  part  r  of  book)  for  -«<x  <  * 

A  forward  difference  equation  is  developed  and  studied  from  the  stand¬ 
points  of  stability,  convergence,  and  discretization  error.  Also  in  part 
2,  for  linear  problems  on  a  finite  interval,  the  forward  difference, 
backward  difference,  Crank-Nicolson,  and  Dufort-Prankel  schemes  are  con¬ 
sidered  with  respect  to  stability,  convergence,  and  discretization  method. 
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35.  Fowler,  C.M.,  "Analysis  of  Numerical  Solutions  of  Transient  Hast  Flow 
Problems,"  Quarterly  of  Applied  Mathematics,  Vol.  3,  p.  361-76,  1946 

Solutions  of  the  difference  equation  for  the  temperatures  Tx  t  for 
one  dimensional  conduction  are  considered.  The  author  uses  contour 
integrals  of  particular  solutions  of  the  difference  equation  to  give 
the  temperatures  in  terms  of  polynomials  and  trigonometric  functions. 

The  convergence  of  his  solutions  to  the  well-known  solutions  of  the 
corresponding  problems  in  partial  differential  equations  as  Ax  and  At* 

0  is  investigated, 

36.  Fox,  L.  (editor),  "Numerical  Solution  of  Ordinary  and  Partial 
Differential  Equations,  "  Addison-Wesley  Publ,  Co.,  Inc.,  Palo  Alto, 
California,  1962,  S09  p. 

Various  finite  difference  schemes  for  partial  differential  equations 
are  discussed  with  regard  to  convergence,  stability  and  computational 
error. 

37.  Fox,  l.,  "Some  Improvements  in  the  Use  of  Relaxation  Methods  for  the 
Solution  of  Ordinary  and  Partial  Differential  Equations,"  Proceedings 
of  the  Royal  Society  of  London,  Series  A,  Vol.  190,  p,  31-S9,  1947 

Boundary  value  problems  associated  with  ordinary  or  partial  differential 
equations  are  commonly  solved  by  the  use  of  difference  equations  which 
are  solved  by  successive  approximations.  Usually  a  derivative  is  replaced 
by  the  leading  term  of  a  finite  difference  series  for  the  derivative  and 
a  small  interval  is  used  to  obtain  the  desired  accuracy.  The  author  pro¬ 
poses  to  use  higher  order  differences  and  a  larger  interval,  which  gives 
a  smaller  number  of  unknowns  to  be  found.  Two  examples  are  worked,  one 
being  Poisson's  equation.  Two  further  examples  show  the  application 
of  the  method  to  curved  boundaries. 

38.  Frankel,  S.P.,  "Convergence  Rates  of  Iterative  Treatments  of  Partial 
Differential  Equations,  Mathematical  Tables  and  Other  Aids  to  Computation, 
Vol.  4,  p.  65-75,  1950 

Convergence  rates  are  estimated  for  several  iterative  methods  of  sol¬ 
ution  for  the  Laplace  and  biharmonic  equations.  The  methods  used  for 
the  Laplace  euqations  are  (1)  Richardson,  (2)  Liebmann.a  ■  1/4, 

(3)  Liebmann,  optimum  a,  (4)  second  order  Richardson.  Quoted  from  the 
author’s  conclusions:  "It  is  thus  seen  that  with  a  fairly  fine  mesh  the 
calculating  time  required  with  the  slower  machines  is  uncomfortably  long 
for  the  Laplace  equation  ...  if  the  normal  Richardson  method  is  used." 

39.  Freed,  N.H.,  and  Rallis,  C.J.,  "Truncation  Error  Estimates  for 
Numerical  and  Analog  Solutions  of  the  Heat  Conduction  Equation," 

ASME  Journal  of  Heat  Transfer,  Vol.  83,  p.  382-3,  1961 

The  authors  describe  a  method  for  obtaining  an  estimate  of  truncation 
error  for  fully  finite  difference  forms  of  the  heat  conduction  equation. 

It  may  be  used  with  manual  and  analog  methods  if  the  error  due  to  mesh 
size  is  relatively  large.  Error  estimates  for  a  case  of  one  dimensional 
flow  are  derived  in  an  example,  using  the  backward  difference  equation. 

40.  Frocht,  M.M.,  "A  New  Approach  to  the  Numerical  Solution  of  Laplace's 
Equations,"  Numerical  Methods  of  Analysis  in  Engineering,  p.  75-80 
The  MacMi Ilian  Co,,  N.Y.,  1949. 
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The  author  discusses  a  procedure  for  obtaining  good  initial  values 
for  interior  grid  points  which  is  applicable  to  iterative  and  relaxation 
methods.  The  procedure  requires  known  boundary  values  and  is  tensed  the 
Linear  Rosette  Method. 

41.  Causer,  G.R.,  "The  Stability  of  Three  Finite  Difference  Methods  of 
Solving  for  Transient  Temperatures,"  presented  at  the  Fifth  U.S.  Navy 
Symposium  on  Aeroballistics  on  18  Oct.  1961,  23  p. 

The  paper  contains  the  results  of  stability  investigations  for  three 
types  of  finite  difference  equations;  which  are  the  forward,  mid  and 
backward  difference  types.  For  each  type,  the  or.e  dimensional  heat 
balance  equations  are  presented  for  three  cor* inat ions  of  heat  transfer 
modes,  which  are  (1)  conduction,  (2)  conduction  and  convection,  (3)  con¬ 
duction,  convection  and  radiation.  A  stability  criterion  is  developed 
for  each  type  of  equation  with  each  combination  of  modes. 

42.  Cill,  J.,  "A  Process  for  the  Step-by-Step  Integration  of  Differential 
Equations  in  an  Automatic  Digital  Computing  Machine,"  Proceedings  of 
the  Cambridge  Philosophical  Society,  Vol,  47,  p,  96-108,  19S1 

The  article  presents  a  modified  4th  order  Runge-Kutta  process  for  systems 
of  1st  order  ordinary  differential  equations.  The  process  described 
requires  a  small  number  of  storage  spaces  for  each  integration  step.  The 
effect  of  truncation  and  round-off  error  is  discussed  and  illustrated  by 
a  numerical  example. 

43.  Goodwin,  E.T.,  Clenshaw,  C.W.,  Martin,  C.h'.,  Miller,  G.F.,  Olver,  F.W.J., 
and  Wilkinson,  J.H.,  "Modern  Computing  Methods,  ",  Philosophical  Library, 
N.Y.,  170  p.,  1961 

The  book  includes  five  chapters  on  matrices,  two  each  on  ordinary 
and  partial  differential  equations,  and  one  on  finite  difference 
methods.  Comments  are  made  on  adaptation  to  desk  and  automatic 
computation  and  careful  attention  is  paid  to  the  assessment  of  comp¬ 
utational  error  for  the  methods  described. 

44.  Graybeal,  G.  E,,  "Time  and  Total  Truncation  Error  Analysis  in  Heat 
Transfer  Calculations",  Master's  Thesis,  West  Virginia  University,  1963, 

59  p. 

A  method  is  presented  for  evaluating  time  and  total  truncation 
errors  encountered  in  a  finite  difference  solution  of  problems  defined 
by  a  parabolic  partial  differential  equation.  The  example  problem  is 
a  one-dimensional  heat  transfer  situation  with  a  step  function  temperature 
change  on  the  surface.  Curves  are  presented  to  aid  in  the  evaluation 
of  time  truncation  and  total  truncation  error. 

45.  Greenuood,  J.  A.,  "Implicit  Numerical  Methods  for  the  Heat  Conduction 
Equation,"  British  Journal  of  Applied  Physics,  Vol.  13,  No.  11,  p.  571-2 

The  author  shows  that  the  Liebnann  finite  difference  approximation  is 
to  be  preferred  to  the  Crank-Nicolson  form  in  certain  cases  involving 
variable  diffusivity. 

46.  Hamming,  R,  W.,  "Numerical  Methods  for  Scientists  and  Engineers," 
McGraw-Hill  Book  Co,.,  Inc.,  San  Francisco,  1962,  411  p. 

The  book  includes  a  discussion  of  difference  calculus  and  difference 
equations.  Functions  of  mo:e  than  one  variable  are  not  treated. 


r 


95 


NWC  TP  5143 


47,  Hawkins,  G.A.,  and  Agnew,  J.  T,,  "Solution  of  Transient  Heat  Conduction 
Problems  by  Finite  Differences,"  Purdue  Univ,  Eng,  Experiment  Station 
Research  Series,  No.  96,  1947,  38p, 

The  purpose  of  the  authors  is  to  bring  the  treatment  for  slabs, 
cylinders,  and  spheres  together.  Analytical  methods  are  discussed  in 
detail  for  unidirectional  flow  in  slabs  and  radial  flow  in  cylinders, 

48,  Henrici,  P.,  "Elements  of  Numerical  Analysis,"  J,  Wiley  and  Sons, 

Inc,,  N.  Y.,  328  p.,  1964. 

This  is  an  introductory  book  and  it  presents  a  small  range  of  subjects 
with  thorough  coverage  on  each  rather  than  treating  a  large  number  of 
techniques  superficially.  It  contains  much  material  on  difference 
equations.  Difference  methods  for  ordinary  differential  equations  are 
developed  using  Taylor's  series.  The  Runge-Kutta  method  is  treated 
only  briefly. 

49.  Herriot,  J.G.,  "Methods  of  Mathematical  Analysis  and  Computation," 

J.  Wiley  and  Sons,  Inc.,  N.Y.,  198  p.  1963 

This  book  is  intanded  for  use  by  engineers,  tnd  it  deals  with  only 
the  best  known  numerical  methods.  The  emphasis  is  on  methods  suitable 
for  use  on  high  spaed  computers.  The  subjects  covered  include:  Hum* 
erical  differentiation  and  integration,  roots  of  equations,  solution 
of  simultaneous  linear  equations,  solution  of  ordinary  differential 
equations,  solution  of  partial  differential  equations. 

50.  Hill,  P.R.,  "A  Method  of  Computing  the  Transient  Temperature  of  Thick 
Walls  from  Arbitrary  Variation  of  Adiabatic  Wall  Temperature  and  Heat 
Transfer  Coefficient,"  NACA  Tech.  Note  410S,  St  p.,  Oct.  19S7 

Quoted  from  the  author's  introduction  "...  simple  method  is  developed 
for  the  calculation  of  the  temperature  history  of  the  surfaces  of  a 
thick  wall  or  of  any  plane  within  the  wall.  The  procedure  is  to  select 
from  a  table  a  set  of  coefficients  which  depend  on  the  physical  properties 
of  the  wall.  These  coefficients  and  other  data  are  substituted  into 
explicit  algebraic  formulas  to  determine-  the  temperature  of  the  heated 
wall  surface.  If  the  heat  transfar  coefficients  are  known,  no  guess 
or  iteration  procedure  is  required.  ....  For  equal  time  step  sizes, 
the  method  is  more  accurst*  than  more  laborious  numerical  methods," 

The  method  uses  concepts  called  'time  series'  and  '  unit  triangle 
variation  of  surface  temperature', 

51.  Hyman,  M.A.,  "On  the  Numerical  Solution  of  Partial  Differential  Equations," 
Thesis,  Technishe  Hogeschool  te  Delft,  108  p,,  1953, 

The  paper  consists  of  four  chapters  and  an  appendix.  Stress  is  placed 
on  methods  suitable  for  use  with  computers.  Ch,  1  treats  convergenca 
and  stability  of  difference  equation  solutions;  Chp.  2  treats  convergence 
and  extrapolation  of  difference  solutions  for  parabolic  equations, 

Ch,  3  discusses  elliptic  equations;  Ch.  4,  is  on  hyperbolic  equations. 

52,  Juncosa,  M.L.,  and  Young,  D.M.,  "On  the  Convergence  of  a  solution  of  a 
Difference  Equation  to  a  Solution  of  the  Equation  of  Diffusion," 

Proceedings  of  the  American  Mathematical  Society,  Vol.  S,  p,  168-74,  19S4. 

Several  sharp  convergence  theorems  are  proved.  The  equation  treated  in 
the  principal  theorem  is  3u/3t  ■  32u/3x2.  The  approximation  used  is: 
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U,  fe+Aftfl.-.M,  *ith 

0<6t/6x2  <1/2,  MAX  ■  1.  Um  (X,t)  is  shown  to  converge  uniformly  to  u(x,t) 
in  the  region  0_<X<1  t>to>0  as  sr*", 

S3.  Juncosa,  M.L.,  and  Young,  D.M.,  "On  the  Order  of  Convergence  of  Solutions 
of  a  Difference  Equation  to  a  Solution  of  the  Diffusion  Equation,"  Journal 
of  the  Society  for  Industrial  and  Applied  Mathenatics,  Vol.  1,  p.  111-35, 
1953 

For  3u/3t  ■  32u/3x2  Is  substituted 


W 1  J£L  r.r.feiti. 

at 


v  (x  ♦  Ax.t)  -  2v(x,t)  •»  vfx-Ax.t). 


The  difference  solution  is  discussed  with  regard  to  attempts  to  improve 
the  solution  by  "extrapolation  to  zero  grid-size." 


54.  Kaplan,  B,  and  Clark,  N.,  "Accuracy  and  Convergence  Techniques  for 
Implicit  Numerical  Solution  of  the  Diffusion  Equation  for  Transient 
Heat  Transfer,"  Transactions  of  the  American  Nuclear  Society,  Vol, 
4,  No,  1,  p.  80-81,  June  1961 


55.  Kardas,  A.,  "Errors  in  a  Finite-Difference  Solution  of  the  Heat  Flow 
Equation,  ASME  Journal  of  Heat  Transfer,  Vol.  86,  p.  561-2,  1964 


Wiis  note  gives  magnitudes  of  discretization  errors  incurred  in  a 
finite  difference  solution  of  the  heat  flow  equation  in  a  symaetric 
slab  with  the  boundary  conditions  of  the  third  kind."  author’s  abstract 

56.  Kunz,  K.S.,  "Numerical  Analysis,"  McGraw-Hill  Book  Co.,  Inc,,  New 
York,  381  p.,  1957. 

The  book  was  written  for  engineers.  It  includes  numerous,  illustrative 
examples.  Iterative  methods  receive  only  cursory  attention.  The  book 
includes  chapters  on  ordinary  and  on  partial  differential  equations. 

An  appendix  on  the  estimation  of  errors  in  numerical  computation  is 
included, 

57.  Lance,  G.N.,  "Numerical  Methods  for  High  Speed  Computers,"  Iliffe  and 
Sons,  Ltd.  London,  166  p.,  1960 

General  text  with  descriptions  of  methods,  often  brief,  some  widely 
used  methods  are  omitted.  Discusses  Runge-Kutta  methods  and  description 
of  Aitkerfs  62  process. 


58.  Lapidus,L "Digital  Computation  for  Chemical  Engineers,"  McGraw-Hill 
Book  Co.,  Inc.,  N.Y,  406  p.,  1962, 


General  text  which  includes  polynomial  approximation,  ordinary  and 
partial  differential  equation,  matrix  solution  of  systems  of  linear 
algebraic  equation,  and  etc,.  Contains  a  description  of  the  Tridag 
method. 


59,  Larkin,  B.K.,  "Some  Finite  Difference  Methods  for  Problems  in  Transient 

Heat  Flow,"  Chemical  Engineering  Progress  Symposium  Series,  Vol,  61,  No.  59 
1965,  p.  1-11 

Four  explicit  methods  for  digital  solution  of  transient  heat  flow  are 
compared.  The  superior  stability  of  the  newer  methods  is  noted.  The 
discussion  is  relevant  to  the  design  of  booster  vehicles,  launch  systems, 
vehicles  for  space  travel,  and  re-entry  heat  shields. 
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60.  Lax,  P.D, ,  and  Richtmyer,  R.D.,  "Survey  of  the  Stability  of  Linear 
Finite  Difference  Equations,"  Communications  on  Pure  and  Applied 
Mathematics,  Vol.  9,  p.  267-93,  19S6 

The  paper  treats  the  numerical  solution  of  initial  value  problems 
by  finite  difference  methods  by  a  Sequence  of  calculations  with 
increasingly  finer  mesh.  The  question  is  whether  the  solution  converges 
to  the  true  solution  of  the  initial  value  problem  as  the  mesh  is  made 
finer.  A  stability  definition  is  given  in  terms  of  the  uniform  bounded¬ 
ness  of  a  certain  set  of  operators,  and  it  is  shown  that,  with  suitable 
circumstances,  for  linear  initial  value  problems,  stability  is  necessary 
and  sufficient  for  convergence  in  a  certain  uniform  sense  for  arbitrary 
initial  data.  Two  different  approximations  to  the  heat  equation  are 
considered,  one  being  a  general  two  level  formula  and  the  other 
the  DuFort-Frankel  equations, 

61.  Lea,  R.N.,  "Stability  of  Multistep  Methods  in  Numerical  Integration," 

NASA  TN  D-2772,  Manned  Spacecraft  Center,  Houston,  Texas,  April,  1965,  16  p. 

The  paper  contains  a  discussion  of  the  stability  of  solutions  of 
differential  equations  obtained  by  using  difference  equations.  An 
original  development  leading  to  a  definition  of  stability  shows  a 
relationship  between  stability  and  certain  properties  of  the  differential 
equation  t»  be  solved.  The  example  worked  in  the  paper  is  an  ordinary 
differential  equation. 

62.  Lees,  M. ,  "Approximate  Solutions  of  Parabolic  Equations,"  Journal  of 
the  Society  for  Industrial  and  Applied  Mathematics,  Vol.  7,  p.167-83,1959 

The  author  discusses  the  convergence  of  numerical  solutions  of  partial 
differential  equations.  The  analysis  is  based  on  energy  methods.  The 
author  derives  a  priori  bounds  for  solutions  of  linear  parabolic  difference 
equations,  then  applies  them  to  establish  the  convergence  of  a  difference 
solution  to  a  non-linear  parabolic  equation,  Crank-Nicolson  type 
difference  equations  are  also  treated. 

63.  Lees,  M.,  "A  Priori  Estimates  for  the  Solutions  of  Difference  Approximations 
to  Parabolic  Partial  Differential  Equations,"  Duke  Mathematical  Journal, 
Vol.  27,  p.  297,  311,  1960. 

The  author  derives,  using  energy  methods  a  priori  estimates  for  the 
solutions  of  several  difference  analogs  of  parabolic  partial  differential 
equations.  All  standard  two  ievel  difference  equations  are  discussed  and 
two  simple  three  level  formula  are  also  treated.  Arguments  are  presented 
in  detail  for  the  heat  equation,  and  generalizations  are  indicated. 

64.  Leppert,  G, ,  "Stable  Numerical  Solution  for  Transient  Heat  Flow,"  ASME 
Paper  No.  53-F-4;  also  published  Amer,  Soc.  Naval  Engrs.  Journal, 

Vol.  6S,  No.  4.  Nov,  1953,  p.  741-52 

An  implicit  finite  difference  formula  for  numerical  integration  of  the 
conduction  equation  is  described.  It  is  shown  to  offer  a  computing 
time  saving  over  previously  used  methods.  The  method  is  a  simple 
algebraic  procedure  for  use  on  a  desk  calculator,  which  removes  the 
necessity  for  iteration  or  substitution  at  each  time  step. 

65.  Leutert,  W.,  "On  the  Convergence  of  Approximate  Solutions  of  the  Heat 
Equation  to  the  Exact  Solution,  "Proceedings  of  the  American  Mathematical 
Society,  Vol.  2,  p.  433-39,  1951 
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The  author  discusses  O’Brien,  Hyman,  and  Kaplan's  criticisms  of  the 
Richardson  difference  equation  (criticises  contained  in  "A  study  of  the 
Numerical  Solution  of  Partial  Differential  Equations"  by  above  named 
authors).  He  shows  that  the  Richardson  difference  equation  is  always 
convergent  if  the  initial  values  of  the  solution  are  chosen  in  a  specifid 
way. 

66.  Leutert,  W, ,  "On  the  Convergence  of  Unstable  Approximate  Solutions  of  the 
Heat  Equation  to  the  Exact  Solution,"  Journal  of  Mathematics  and 
Physics,  Vol.  30,  p.  245-SI,  1952. 

The  differential  equation  considered  is  du/dt  ■  du2/dx2.  The  finite 
difference  equation  considered  is  v(X,t+dt)  -  v(X,t)  ■  r[v(x*Ax,t)  ♦ 
\(X-AX,t)  -2v(x,t)],  where  (AX)2r  »A  t.  It  is  known  that  for  r»l/2  the 
solution  of  the  finite  difference  equation  is  unstable,  The  author 
proves  that,  even  so,  there  exist  for  every  fixed  r,  solutions  of  the 
difference  equation  that  converge  to  the  solution  of  the  differential 
problem  as  AX-»0. 

67.  Liebmann,  G.,  "Solution  of  Transient  Hest  Flow  and  Heat  Transfer 
Problems  by  Relaxation,"  British  Journal  of  Applied  Physics,  Vol,  6, 

No.  4,  Apr.  19SS,  p.  129-35 

This  illustrates  that  by  choosing  a  suitable  finite  difference 
approximation,  parabolic  partial  differential  equations  can  be  converted 
into  a  series  of  boundary  value  problems  of  Poisson  type,  which  can  be 
solved  by  the  Southwell  relaxation  technique.  A  very  stable  solution 
is  obtained  for  all  values  of  the  time  interval  by  using  a  backward 
difference  approximation. 

68.  Lotkin,  M.,  "On  the  Accuracy  of  Runge-Kutta's  Method",  Mathematical 
Tables,  and  Other  Aids  to  Computation,  Vol.  5,  p.  128-33,  1951 

The  author  obtains  a  bound  for  the  error  in  Kutta’s  fourth  order 
method  (generally  know  as  the  Runge-Kutta  method). 

69.  Lotkin,  M.,  "The  Numerical  Integration  of  Heat  Conduction  Equations," 
Journal  of  Mathematics  and  Physics,  Vol.  37,  p.  178-87,  1958 

The  author  discusses  difference  equation  approximations  to  the 
equations  of  unsteady  one-dimensional  heat  conduction  in  composite 
media:  ♦(!<>  3u/3f  »3/3x  (kl(*»)  u/3“.)  where  k^  and  $(»)  denote 
known  functions  of  u.  The  convergiflce  is  established  for  kC*0  ■  constant 
and  the  convergence  rate  is  estimated.  An  example  is  given,  and  numerical 
data  compares  the  approximate  and  exact  solutions. 

70.  Luke,  Y.L.,  "Numerical  Solution  of  the  Heat  Conduction  Equation," 

Chemical  Engineering,  Vol.  68,  No.  1,  Jan.  9,  1961,  p.  95,  102 

The  article  discusses  the  numerical  integration  of  the  heat  conduction 
equation,  and  the  computation  of  flux  and  heat  transfer.  Included  are 
tables  listing  numerical  solutions  for  constant  thermal  coefficients. 
Variable  diffusivity  is  also  discussed. 

71.  Lynn,  L.L.,  and  Meyer,  J.E.,  "A  Numerical  Comparison  of  the  Implicit 
and  Explicit  Techniques  for  the  Convective  Boundary  Condition,"  ASME 
Journal  of  Heat  Transfer,  Vol.  85,  p,  280 r  81,  1963, 

The  article  compares  the  results  of  the  Crank-Nicolson  implicit  method 
with  those  of  Back's  explicit  method  for  1000°F  step  change  in  ambient 
temperature.  Some  comments  and  conclusions  are:  (1)  The  implicit 
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calculations  for  surface  temperature  are,  in  most  cases,  sore  accurate 
than  the  Bach  method  for  a  given  C1m  step  site,  (2)  The  implicit 
method  permits  variable  sesh  spacing  with  no  concern  for  stability. 

(3)  The  use  of  Gauss  elimination  in  the  implicit  method  involve*  more 
effort  than  the  explicit  method  per  point  per  time  step,  but  a  dacreaae 
in  the  nuaber  of  time  steps  needed  for  a  given  accuracy  often  offsets  this, 

72.  Hacen,  N« ,  "Numerical  Analysis,"  John  Wiley  and  Sons,  Inc.,  K.Y.,  1963, 

161  p. 

In  this  book,  the  emphasis  is  on  methods  for  use  with  digital  computers. 

Some  of  the  topics  covered  are:  Iterative  methods  for  solving  equations; 
matrices  and  systems  of  linear  equations;  difference  equations. 

73.  Mann,  W.R.,  Bradshaw,  C.L. ,  Cox,  J.G.,  "Improved  Approximations  to 
Differential  Equations  by  Difference  Equations,”  Journal  of  Mathematics  and 
Physics,  Voi.  3S,  p.  408-15,  1957 

This  article  shows  that  the  truncation  error  in  using  a  difference 
equation  to  approximate  a  differential  equation  can  be  reduced  by 
modifying  the  coefficients  of  the  difference  equation  from  those  normally 
used.  The  difference  equation  is  expanded  In  a  Taylor's  series  and 
the  expansion  compared  with  the  leading  terms  of  the  original  difference 
equation,  the  comparison  provides  a  correction  which  allows  a  reduction 
in  truncation  error  without  increasing  the  order  of  the  difference 
equation. 

74.  Martin,  D.K.,  "Runge-Kuttt  Methods  for  Integrating  Differential  Equations 
on  High  Speed  Digital  Computers,”  The  Computer  Journal,  Voi.  1,  p,  118*123 
1958 

The  author  describes  three  adaptations  of  Runge-Kutta  procedures  for 
ordinary  differential  equations,  due  to  Gill,  Strachey,  and  Boulton; 
and  he  proposes  an  alternate  method  devised  to  save  storage  space  and 
based  on  Xutta's  Simpson  rule  method.  Comparative  errors  and  computational 
experience  with  the  various  methods  are  described. 

?S.  Milne,  W.E„,  "Note  on  the  Runga-Kutta  Method,”  Journal  of  Research, 

National  Bureau  of  Standards,  'J44,  p.  549*50,  1950 

Comparison  of  the  Runge-Kutta  method  with  a  step-by-step  tethod  of 
numerical  quadrature  shews  the  Runge-Kutta  technique  to  be  much  less 
accurate  in  seme  cases. 

76.  Milne,  W.E. ,  "Numerical  Solution  of  Differential  Equations,”  J,  Wiley  and 
Sons,  Inc.,  New  York,  27S  p,,  1953, 

The  book  offers  general  coverage,  Including  both  ordinary  and  partial 
differential  equations.  Notes  on  large  scale  computers  occur  in  an 
appendix, 

77.  Milne,  W.E.,  "Numerical  Methods  Associated  With  Laplace's  Equation," 
Proceedings  of  a  Second  Symposium  on  Large  Scale  Digital  Calculating 
Machinery,  1949,  p.  152-63,  Harvard  University  Press,  Cambridge,  Mass., 1951 

This  is  a  review  of  some  difficulties  which  occur  in  solving  partial 
differential  equations  b>  the  method  of  differences,  using  large  scale 
digital  machines. 

78.  Mitchell,  A.P,,,  "Round-Off  Errors  in  Solution  of  Heat  Conduction 
Equation  by  Relaxation  Methods,"  Applied  Science  Research,  Sect.  A, 

Vol,  4  n.  2,  1953,  p.  109-19 


100 


NWC  TP  5143 


A  method  is  developed  for  assessing  the  magnitude  of  round-off  errors, 
a  stable  six  point  finite  difference  approximation  is  used,  which  is 
relaxationsl  in  distance  coordinate  and  step  by  step  in  tine  coordinate 
of  the  residuals,  Formulas  are  derived  for  obtaining  round-off  errors 
for  several  different  distributions  of  residuals. 

79.  Mitchell,  A.R.,  Fairweather,  G.,  "Improved  Fores  of  the  Alternating 
Direction  Methods  of  Douglas,  Peacenan,  and  Rachford  for  solving 
Parabolic  and  Elliptic  Equations,"  Nuaerische  Mathematik,  Vol,  6, 
p.  285-92,  1964 

For  the  heat  conduction  equation  Ut  ■  Uxx  ♦  U„y  on  a  rectangular 
domain,  the  authors  derive  generalizations  of'  the  Peacenan-Rachford 
alternating  direction  difference  schemes.  A  special  choice  of  para¬ 
meters  leads  to  a  stable  scheme  with  fourth  order  accuracy, 

80.  Mitchell,  D.3.,  "An  Error  Analysis  of  Numerical  Solutions  of  the 
Transient  Heat  Conduction  Equation,"  Master's  Thesis,  Rept.  No. 

GA/PU/6S-10,  AF  Inst,  of  Tech,,  Wright-Patterson  AFB,  Ohio,  Aug,  1965, 

111  p.  AD-621  274 

The  paper  presents  a  comparison  of  the  Crank-Nicolson  and  Crandall 
methods  in  finding  transient  temperatures  in  a  semi-infinite  slab 
with  convection.  The  results  are  compared  with  an  exact  analysis.  The 
Crandall  method  gives  more  accurate  results  than  the  Crank-Nicolson 
under  close  node  spacing  conditions.  Accuracy  improvement  factors 
are  determined  for  the  two  methods. 

81.  Muchnik,  G,  F.,  "Solution  of  Heat  Conduction  Problems  by  the  Grid 
Method,"  NASA  TTF-151,  April  1964,  15  p,  (translated  from  Russian). 

An  approximate  method  is  suggested  for  solving  heat  conduction  problems 
using  the  grid  method.  The  temperature  of  any  point  is  related  to  the 
temperatures  of  other  points  by  coefficients  of  relationship  or 
"weights",  which  do  not  depend  or.  the  boundary  conditions.  The  weights 
are  found  by  a  finite  difference  method.  The  author  claims  this  method 
to  be  simpler  and  more  exact  than  the  usual  finite  difference  method. 

82.  Murray,  W.D.,  and  Landis,  Fred,  "The  Effect  of  Spacewise  Lumping  on 
the  Solution  Accuracy  of  the  One-Dimensional  Diffusion  Equation,"  ASME 
Journal  of  Applied  Mechanics,  Vol.  29,  p,  629-36,  1962, 

The  authors  evaluate  the  truncation  errors  inherent  in  a  spacewise 
difference  formulation  of  the  one-dimensional  heat  diffusion  equation 
under  general  boundary  conditions.  The  error  between  the  semi- 
dlacrete  end  exact  solutions  is  evolved  by  matrix  algebra  and  the 
Laplace  transform.  An  illustration  shows  the  errors  for  the  case  of 
a  syMtatrically  heated  slab, 

83.  Nielsen,  K.L.,  "Methods  in  Numerical  Analysis,"  The  MacMillan  Co., 

New  York,  382  p,,  19S6. 

The  author  emphasizes  methods  suitable  for  desk  calculating. 

One  chapter  treats  ordinary  and  partial  differantial  equations,  but  it 
presents  only  an  outline  of  the  topic.  The  methods  of  Euler,  Milne, 
and  Runge-Kutta  are  given  for  ordinary  differential  aquations,  Liebmann's 
method,  relaxation,  and  step  by  step  methods  are  given  for  partial  differen¬ 
tial  equations. 

84.  O'Brien,  G.G.,  Hyman,  M.A, ,  and  Kaplan,  5.,  "A  Study  of  the  Numerical 
Solution  of  Fartial  Differential  Equations,"  Journal  of  Mathematics 
and  Physics,  Vol.  29,  p.  223-51,  1951 
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This  article  is  a  discussion  of  methods  for  the  analysis  and  improvement 
of  the  stability  of  finite  difference  equations  used  in  the  numerical 
solution  of  partial  differential  equations.  A  discussion  of  convergence 
is  included.  The  heat  conduction  equation  of  the  form3+/3t  «32+/3x‘ 
is  used  as  an  illustration.  The  equivalent  difference  forms  of  Richardson, 
von  Neumann  and  Hartree  are  discussed. 

85.  Peaceman,  D.W.,  and  Rachford,  H.H.,  "The  Numerical  Solution  of  Parabolic 
and  Elliptic  Differential  Equations,"  Journal  of  the  Society  for 
Industrial  and  Applied  Mathematics,  Vol,  3,  p.  28-41,  195S. 

The  authors  treat  Ut  all**  ♦  “yy.  by  implicit  difference  mothods.  They 
develop  a  process  which  is  staoie  for  all  mesh  ratios  ((Ax)2/t).  They 
discuss  the  solution  over  a  square,  showing  the  work  saving  for  their 
method  compared  to  usual  ones,  For  elliptic  equations  their  method  is 
a  form  of  line  relaxation  rather  than  the  usual  point  relaxation  schemes. 

86.  Plunkett  R.,  "On  the  Rate  of  Convergence  of  Relaxation  Methods," 

Quarterly  of  Applied  Mathematics,  V.  19,  p.  263-66,  1952. 

This  article  compares  the  convergence  rate  for  the  relaxation  method 
of  solving  partial  differential  equations  to  the  results  obtained  by 
Frankel  (Math.  Tables  and  Other  Aids  tc  Comp.  v.  4,  p.  65-75,  1950)  for 
an  iteration  method.  It  is  concluded  that  the  relaxation  method  gives 
no  saving  for  the  Poisson  and  biharmonic  equations  and  is  more  difficult 
to  program  then  the  iterative  method, 

87.  Poppe,  R.T.,  "An  Investigation  of  Convergence  Techniques  for  Implicit 
Numerical  Solution  of  the  Diffusion  Equation  for  Transient  Heat  Transfer," 
Master's  thesis,  AFIT/GA/phys/63-8,  AF  Inst,  of  Technology,  Wright-Patterson 
AF  Base,  Ohio,  Aug.  1963,  163  p.,  AD-419  310. 

This  paper  contains  the  results  of  an  investigation  of  two  techniques 
for  increasing  the  rate  of  convergence  of  the  Causs-Siedel  method  of 
implicit  numerical  solution  of  the  diffusion  equation  of  transient  heat 
flow.  A  sample  problem  is  solved  to  provide  the  necessary  comparison. 

The  results  provide  a  theoretical  basis  for  the  adapted  Wegstein  technique. 
This  theoretical  basis  brings  to  light  the  fact  that  successive  over¬ 
relaxation  and  the  adapted  Negstein  technique  are  based  on  same  theoretical 
background.  A  procedure  based  on  estimating  the  maximum  eigenvalue  of 
the  method  of  successive  displacements  is  used  to  make  an  approximation 
of  the  relaxation  factor  for  successive  over- relaxation. 

88.  Prager,  W. ,  "Introduction  to  Basic  FORTRAN  Programming  and  Numerical 
Methods"  Blaisdell  Publishing  Company,  New  York,  202  p.,  1965 

This  book  contains  information  on  Aitkcn's  62method  and  Steffensen's 
method  for  accelerating  the  convergence  of  iterative  processes. 

89.  Price,  P.H.,  Slack,  M.R.,  "Stability  and  Accuracy  of  Numerical  Solutions 

of  Heat  Plow  Equation,"  British  Journal  of  Applied  Physics,  Vol.  3,  No.  12 
Dec.  1952,  p.  379-84 

The  authors  describe  a  new  method  of  deriving  stability  conditions. 

Its  application  to  heat  conduction  with  variable  thermal  diffusivity 
and  heat  transfer  by  convection  at  surface  is  treated.  A  new  finite 
difference  representation  of  surface  heat  flux  equation  is  given. 

90.  Price,  P.H, ,  and  Slack,  M.R, ,  "Effect  of  Latent  Heat  on  Numerical 
Solutions  of  Heat  Flow  Equation,"  British  Journal  of  Applied  Physics, 

Vol.  5,  No.  8,  Aug.  1954,  p.  285-7 
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The  article  treats  the  stability  and  accuracy  of  finite  difference 
solutions  of  the  heat  flow  equation  with  latent  heat  evolution.  A 
dimensionless  group  is  developed  which  governs  the  -appearance  of  in¬ 
accuracies  peculiar  to  numerical  solutions  involving  latent  heat. 

91.  Ralston,  A.,  Wilf,  H.S.,  "Mathematical  Methods  of  Digital  Computer," 

John  Wiley  £  Sons,  Inc.,  Mew  York,  1960 

A  section  on  the  solution  of  linear  equations  includes  discussions  of 
Causs-Seidel  iteration,  the  conjugate  gradient  method.  Matrix  Inversion 
by  Rank  Annihilation,  Matrix  Inversion  by  Monte  Carlo  Methods  -  flow 
chart  and  formula  for  estimating  running  time  arc  included  for  each  method. 

92.  Round,  O.P.,  Newton,  R.,  and  Redberger,  P.J.,  "Variable  Mesh  Si2e  in 
Iteration  Methods  of  Solving  Partial  Differential  Equations  and 
Application  to  Heat  Transfer,"  Chemical  Engineering  Progress 
Symposium  Series,  Vol.  58,  No.  37,  1962,  p.  29-42 

The  article  contains  descriptions  of  some  variable  mesh  systems. 
Computations  are  carried  out  for  steady  state  heat  transfer  from  a 
buried  cylindrical  heat  source.  The  same  accuracy  is  achieved  in  a  shorter 
time  than  with  a  square  mesh. 

93.  Saul'ev,  V.K.,  "Integration  of  Parabolic  Equations  by  the  Method  of 
Nets,"  translated  frost  Russian; MacMillan,  N.Y.,  1964,  Russian  publication 
date  1960. 

There  are  two  sections,  the  first  dealing  primarily  with  stability  and 
convergence  and  the  second  dealing  with  implicit  methods  and  techniques 
for  solving  the  algebraic  equations.  In  the  first  section,  a  number  of 
mesh  schemes  are  considered  and  compared.  The  second  section  gives  a 
fairly  complete  compendium  of  methods  for  solving  linear  algebraic  systems. 

94.  Scarborough, J ^'Numerical  Mathematical  Analysis,"  5th  Edition,  The  Johns 
Hopkins  Press,  Baltimore,  594,  p,,  1962 

There  is  a  chapter  on  partial  differential  equations  which  treats 
difference  quotients,  difference  equations,  the  solution  of  differential 
equations  by  iteration,  the  inherent  error  in  the  solution  of  difference 
equations,  and  relaxation  methods, 

95.  Schenck^H, Jr., "Fortran  Methods  in  Heat  Flow,"  The  Ronald  Press  eo,, 

N.Y.,  289  p,,  1963 

As  indicated  by  the  title,  the  emphasis  is  on  Fortran  methods,  with 
little  discussion  of  theory.  The  chapter  on  one»di*ensional  transient 
flow  follows  the  method  of  Dusinberre,  but  includes  a  brief  discussion 
of  the  relative  merit  of  the  Liebmann  implicit  method.  Sample  Fortran 
programs  are  presented  (using  the  Dusinberre  method  only.)  The  last 
chapter  of  the  book  is  a  short  discussion  of  accuracy  and  of  solution 
speed. 

96.  Schneider,  P.J.,  "Conduction  Heat  Transfer,"  Addison-Wesley  Publishing 
Co.,  394,  p.  1955 

Chapter  12  of  this  book  treats  the  transient  numerical  solution  of 
conduction  problems  using  the  method  of  Dusinberre.  A  discussion  of 
stability  and  convergence  is  included. 

97.  Schuh,H"Finite  Difference  Method  for  Calculating  Transient  Temperature 
Distributions  Due  to  One-Dimensional  Heat  Flow  in  Simple  and  Composite 
Bodies,  "RAE-Lib/Trans-750,  translated  from  VDI  -  Forschungshoft, 

No.  459,  43.  p. 
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The  finite  difference  method  is  Modified  for  finding  transient  temperatures 
due  to  heat  flow  normal  or  parallel  to  the  surface  of  thick  plates  and 
thin  walled  bodies  of  high  conductance.  Details  of  the  Method  are 
discussed  extensively.  The  use  of  high  speed  coaputers  for  solving 
finite  difference  equations  is  discussed. 

98.  Southwell,  R.V.,  "The  Quest  for  Accuracy  in  Coaputations  Using  Finite 
Differences,"  Nuaerical  Methods  of  Analysis  in  Engineering,  p.66-74, 

The  MacMillan  Company,  New  York,  1949. 

The  author  concludes  that  the  bast  accuracy  is  obtained  by  reducing 
the  interval  size  rather  than  using  higher  order  difference  equations 
and  illustrates  a  labor  saving  device  which  facilitates  "advance  to  a 
finer  net". 

99.  Southworth,  R.W.,  DeLeeuw,  S.L.,  "Digital  Coaputation  and  Numerical 
Methods,"  McGraw-Hill  Book  Co.,  N.Y.,  S08  p,,  1965. 

This  book  is  intended  for  use  as  a  textbook  in  a  course  conbining 
FORTRAN  programming,  nuaerical  methods,  and  engineering  applications. 

Chapters  2,  5,  and  4  deal  with  programing,  and  the  reaainder  of  book 
is  concerned  with  nuaerical  methods.  Problems  are  given  at  the  end  of 
each  chapter,  and  engineering  applications  appear  throughout  the  text. 

100.  Stanton,R "Numerical  Methods  for  Science  and  Engineering,"  Prentice-Hall, 
Inc.,  Englewood  Cliffs,  N.J.,  266  p.,  1961 

This  book  is  intended  as  an  undergraduate  introduction  to  numerical 
analysis,  and  is  short  on  precise  theory.  It  stresses  desk  calculator 
methods.  Included  are  discussions  of  ordinary  finite  differences, 
divided  differences,  and  central  differences.  There  is  one  chapter 
on  the  solution  of  differential  equations  by  difference  equation  methods. 

101.  Strang,  W.G.,  "On  the  Order  of  Convergence  of  the  Crank-Nicolson 
Procedure,"  Journal  of  Mathematics  and  Physics,  Vol.  38  p.  141-44, 

1959-60 

The  author  discusses  the  Crank-Nicolson  difference  equation  for  Ut* 

Uxx  ♦  d(x,t).  If  the  solution  of  the  first  boundary  problem  is 
sufficiently  smooth,  the  solution  of  the  difference  equation  converges 
point-wise  with  error  0(  (Ax)2),  if  At  •  0  (AX).  The  proof  makes  use 
of  explicitly  known  eigenfunctions  of  the  process, 

102.  Thomas,  L.H.,  "Numerical  Solution  of  Partial  Differential  Equations 
of  Parabolic  Type,"  Proceedings,  Seminar  on  Scientific  Computation, 

Nov.,  1949,  p.  71-78,  International  Business  Machines  Corp.,  New  York, 

N.Y,  1950 

The  article  contains  an  expository  treatment  of  some  problems  in 
the  numerical  solution  of  parabolic  partial  differential  equations  by 
finite  differences.  Thorr  are  three  major  topics:  (A)  stability  of 
the  finite  difference  representation,  (b)  truncation  errors,  (c)  round¬ 
off  errors.  Methods  for  improving  the  stability  and  reducting  truncation 
errors  are  illustrated, 

103.  Thomas,  L.  H.,  "Stability  of  Solution  of  Partial  Differential  Equations," 
Rept,  No.  NOLR-1132,  Naval  Ord,  Laboratory,  White  Oak,  Md,,  1950,  p.83-94 
title  of  report  is  "Symposium  on  Theoretical  Compressible  Flow,  28  June  1949" 
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The  report  is  a  survey  of  the  present  status  of  the  art  of  stability 
analysis  for  finite  difference  equivalents  of  differential  equations 
(both  ordinary  and  partial).  The  author  defines  two  types  of  instability, 
short  range  and  long  range. 

104.  Traub,  J.F.,  "Iterative  Methods  for  the  Solution  of  Equations," 

Prentice-Hall,  Inc.,  Englewood  Cliffs,  N.J.,  310  p.,  1264. 

A  large  nunber  of  iteration  functions  are  described  and  classified 
according  to  the  efficiency  of  the  algorithm  and  the  amount  of  computa¬ 
tional  labor  involved.  Stress  is  placed  on  methods  for  constructing 
iteration  functions  and  on  determining  their  chief  properties.  Much 
attention  is  focused  on  the  computational  aspects  of  the  topic. 

105.  Turner,  L.R. ,  "Improvement  in  the  Convergence  of  Methods  of  Successive 
Approximation,"  Proceedings,  Computation  Seminar,  Dec.  1949,  p.  135-57 
International  Business  Machines,  Corp.,  New  York,  N.Y. ,  1951 

An  exposition  on  the  well  known  procedure  for  improving  the  convergence 
of  an  iteration  procedure  when  the  steps  form  a  geometric  progression. 

106.  Turton,  F.J.,  "The  Errors  in  the  Numerical  Solution  of  Differential 
Equations,"  The  Philosophical  Magazine,  Vol.  28,  p.  359-63,  1939 

The  article  contains  a  detailed  analysis  of  the  errors  caused  by 
(1)  uncertainty  of  initial  values,  (2)  intrinsic  errors  in  formulae 
used  in  the  step  by  step  method,  (3)  round-off  errors  in  (2),  (4) 
random  errors.  The  author's  conclusion  is  that  "to  insure  no  errors 
to  the  desired  number  of  significant  figures....  requires  that  at 
least  two  formulae  be  used,  in  which  the  intrinsic  errors  are  sub¬ 
stantially  different,  to  check  each  other." 

107.  Varga,  R.S.,  "A  Comparison  of  the  Successive  Overrelaxation  Method  and 
Semi-Iterative  Methods  Using  Chebyshev  Polynomials,"  Journal  of  the 
Society  for  Industrial  and  Applied  Mathematics,  Vol.  5,  p.  39-46,  1957 

The  author  shows  that  the  successive  overreiaxation  method  converges 

at  least  as  fast  as  any  somi-iterative  method  associated  with  the 

Jacobi  method,  Causs-Siedel  method,  or  with  the  successive  overreiaxation 

method  itself.  Successive  overreiaxation  requires  only  the  latest 

iterate  at  any  stage,  whereas  3emi-iterative  methods  require  the  simultaneous 

storage  of  several  iterates;  therefore,  the  author  sees  some  advantage 

in  using  successive  overreiaxation,  instead  of  semi-iterativc  methods, 

with  high  speed  computers. 

108.  Wasow,  ><.,  "On  the  Accuracy  of  Implicit  Difference  Approximations  to 
the  Equation  of  Heat  Flow,"  Tech,  Summary  Rept.  No.  MRC-TSR-2, 

Contract  DA11  022  0RD  2059,  Math,  Research  Center,  Univ.  of  Wisconsin, 
Madison,  Wise.,  15  Apr.  1957,  22p.  PB-167  605 

The  author  discusses  the  convergence,  stability  and  truncation  error  of 
implicit  difference  approximations  to  the  initial  value  problem  defined  by 
the  heat  flow  equation. 

109.  Wegstein,  J.H.,  "Accelerating  Convergence  of  Iterative  Methods:" 
Communications  of  Computing  Machinery,  Vol.  1,  No.  6,  p.  9,  1958, 

The  article  describes  a  method  very  similar  to  Aitken's  62  method. 

It  is  emphasized  that  the  method  can  cause  convergence  in  normally 
divergent  cases.  Several  numerical  examples  arc  included. 
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110.  Wilkes,  J.O.,  "Chemical  Engineering  Workshop:  Part  11,  Numerical 
Methods  for  Partial  Differential  Equations,"  paper  presented  at 
ASEE  annual  meeting,  Pullman,  Wash.,  June  1066 

The  author  describes  3  types  of  finite  difference  approximations  for 
partial  differential  equations,  called  the  forward,  backward  and  central 
difference  types.  The  main  part  of  the  paper  is  a  discussion  of  the 
solution  of  a  heat  transfer  problem  using  finite  differences.  The 
FORTRAN  II  program  employed  is  included  in  full,  as  is  a  special 
subroutine  (Tridng)  for  solving  the  system  of  linear  equations  resulting 
from  application  of  the  finite  difference  approximation. 

111.  Zonneveld,  J.A.,  "Automatic  Numerical  Integration,"  Mathematical 
Centre  Tracts,  No.  8,  Mathematisch  Centrum.  Amsterdam,  1964,  110  p. 

The  author  constructs  a  set  of  Runge-Kutta  formulas  suitable  for  automatic 
adjustment  of  step  size,  The  function  being  integrated  is  evaluated  for 
step  size  one  or  two  additional  times,  and  the  size  of  the  step  is 
adjusted  according  to  this  evaluation,  during  the  integration  process. 

Also  included  are  six  ALGOL  60  program;  for  the  integration  of  first  and 
second  order  differential  equations,  and  five  numerical  examples. 
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119 


DERIVATION  OF  THE  LAPLACE  TRANSFORM  SOLUTION 


This  appendix  is  an  application  of  the  Laplace  transform  method 
to  the  solution  of  the  one-dimensional  Fourier  conduction  equation  for 
a  semi-infinite  solid  and  for  an  insulated  slab,  each  with  a  ramp- 
function  boundary  condition. 


SEMI-INFINITE  SOLID 


The  one-dimensional  Fourier  conduction  equation  is 


=  a  ^  at  T  =  0,  t  =  0 

3T  ax2 


The  Laplace  transform  is 


d2T 

a  =  sT  _  C(0) 
dX 


and  for  the  conditions  of  Eq.  1,  it  becomes 

d^T  _  sT  n 

ox2'" 

and  has  the  general  solution, 

T(X,  s)  =  A  exp Ws / a  X  +  B  exp(-  •Js/a.  X 

Since  the  body  is  a  semi-infinite  solid,  t  -*■  0  as  X  -*•  00 
implying  that 

A  =  0 

and  the  solution  reduces  to 

T(X,  s)  =  B  exp(-  \/s/a  X) 


The  coefficient  B  must  be  determined  from  the  surface  conditions, 
At  X  =  0,  tf  =  f(T)  as  shown  in  the  following  figure. 
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*lCf  exp(-  %/s/a  X) 

^ ‘  1  s2  y/sja  +  h/k 


exg.1-  JUk-lL-  exp(-T  s) 
s2(  y/s/a  +  h/k) 


The  inverse  Laplace  transform  of 


exp(-  s/n  X) 
s  (  s/ot  +  h/k) 


k  c 

—  erfc 
h 


fe)-b4xa(tf  t]ciec(w+^t 


The  inverse  Laplace  transform  of 


exp(-\/s7a  X)  !■ 
s  ( \/~i T7a"  +  h/k)  s 


1  *  £  erfc 


(  zkj-  I  “»[  E  x  +  °(tf  T]er£c(^+  Hi (11) 


and  becomes 


h/.o 

-  I. /  sxp[e  H  X]  er£c(275T  +  t'61)  a) 


The  inverse  Laplace  transform  of 

exp(^X)  eXp(~TlS) 
s (Vs / a  +  h/k)  s 


|/  "‘■[iAPyl*1 

*7=0  m  1 

-  {jjJ^  e*P  X  +  (X-T)  erfc  dxj 

•  10(1-12)]  (13) 
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Let 


6  = 


4aA 


/•X2/4aT  /  x23-2\ 

L  er£cWff>  v 


and  the  first  term  of  Eq.  12  becomes 

erfcu/g)  {-  ~7 — )df$ 

which  equals 


..2  fXl/aT 

~  I  -g  ^(l-erfN/H)dg 

J  g=“ 


Integrating  the  preceeding  integral  by  parts  yields 
T  efrcf^r^l-  XVfT'aTr 


XVT7S"  ex|>(-  4S  ) 


+  7?  5?  (  2  B"1/2  “P(-8)<iB 

40 


which,  by  letting  g  =  a  ,  becomes 
T  erfcj 

and  the  first  term  of  Eq.  12  becomes 


Mr  A'  *L/ 2vrTexp('a2)da 


T  + 


jr 

2a 


) erfc(rs)  -x'ff7S"  exp(-  ) 


(14) 


The  second  term  of  Ecj.  12  becomes 
*T  f  a  \2 


/=0  exp[x  +  °§) x  ]  crfc(2vir + 

(t  x)/=0  exp  “(£)  x]  J1  -  erf(dnr + 1 


“  exp 


d> 
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Integrating  by  parts  yields 


«  (c)  axp(t  x)  | axp[“00 

;xp  -I  € 


T  erfcl 


,  hVS 
+  kV*  6XP 


X 

4^  exp 


(2  Jt  +  k  '®) 

(-  4£ 

(-  s x)  *£, r3/2  exp  (■  4) 


A_1/2  exp! 


(15) 


Integrate  the  first  integral  in  Eq,  15  by  parts  and 

,-1/2 


L  X'in  e*P(~  4)  di  ■  m  axpf  fcr ) 

£0  r3/2  e4  4) « 


r 

2a 


Substituting  in  Eq.  15  gives 


,2VSf 


,  2h,/oT 

+  r  V7'  axp 


[c  *  *  i!;) T]  *'c{ 


"(2^) 


(16) 


Let 


Z  = 


2\JaX 

and  Eq.  16  becomes 


a  00  {  eXP  t  X  +  “(k)  T  erfC(2^T  +  k 

2 


,  2h  r~ 

+  V?r  ^  exp 


2VaT 


-(r+1)  erfc(w)l  (17) 


Summing  Eq .  14  and  17  and  multiplying  by  htf/k.2  determines  the 
temperature  distribution  for  0  <  T  <_T],,  which  is 
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t 


-«(£)  exp[fe x  +  “(k)  TJ  erfc(v§f +  1^) 

^  +  exp[_(^Br)  ]  ]■ 

Similarly,  the  temperature  distribution  for  T  >  is 


'  -  fj[[T  +  £ +  S + :  (s)2]  erfc(w) 

-  a  (I)  e,!p[l  X  +  “(I)  T]  erfc(w  +  I  'fiT) 

-V?  (-fHMrfij 


U" 

-3 

1 

-3 

\  +  i. ™  + 1  (k' 

k2] 

)  erfc 

X 

T1  (U 

}  2a  ’  ah  a  \h, 

Va(T-Tx) 

(18) 


i(t)  «l>[£  X  +  «  (j[)  (t-Tj)]  .rfc[jJ5!pT5+ 


J  (x  +  ir)  6xp(-[wptit]  )1 


(19) 


INSULATED  SLAB 

Again  starting  with  the  one-dimensional  Fourier  conduction  equation, 
one  sees  that 

2 

It  =  a  fx? at  t  =  °»T=°  (20) 

and  the  Laplace  transform  is 


125 


N’WC  TP  5143 


-  I  T(X,s)  =  0  (21) 

dX 

A  general  solution  is: 

T(X,s)  =  A  cosh  y/s/a  X  +  B  sinh  >fs/ a  X  (22) 

The  temperature  of  the  bounding  fluid  expressed  in  terms  of  time,  t> 
is 


f(T)  = 


V 


‘1L 


Tu(T) 


(23) 


A  heat  balance  at  the  surface  of  the  slab  yields 

ht 


3t(0,T)  +  ht  (0  ,T)  =  ^f 

ri 


sx  +  ‘  k?:[T,!(T)  •  (T-Ti)u(T-Ti) 


(24) 


The  Laplace  transform  of  Eq.  24  is 

.M  +  ht(0,s)^ 


1  “p(-v)i 

2  "  2 


(25) 


but 


bvSTT,  t(o,s) 


=  A 


Thus  Eq .  25  becomes 


r  n~  +  a  ^  ^f  f  ^ 

-B  vs/a  + 


expJ-Tjs) 


(26) 


Since  the  slab  is  insulated  at  X  =  L,  —  =  0,  and  from  Eq.  22, 

aX 


\/s7a  | A  sinh  <J s fa  L  +  B  cosh  t/s/a  lj  =  0 
Solving  Eq .  26  and  27  simultaneously  yields 


(27) 


A  = 


!!lt| 

!  cosh  yjs/o.  L  \  | 

[i  exH 

;-v)i 

kTj  ! 

^i/s/a  sinh  ’Js/a  L  +  h/kv/cosh  \/s]a  L J 

2  2 

LS  S  i 
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B  *  - 


/ 

kTi  V 


sinh  \/s/ 


kT,  \  Vs/a  sinh  Vs/a  L  +  h/k  cosh  Vs /a  LJ 


i  espKs}| 


Therefore, 
T(X,  s) 


( cosh  Vs/a  L  cosh  Vs /a  X  -  sinh  Vs/a  L  sinh  Vs/ a  X 

=  kT, 


jf  /  c  _ 

rx  ^  VsTcT  sinh  y/s'fa  L  +  h/k  cosh  y/sTa  L 

r  i  exp(-Tis)i 


2  2 
s  s 


(28) 


which  can  be  put  in  the  following  form: 


T(X,  s)  = 

*1 


cosh  Vs /a  (L  -  X) 


k/h  Vs / a  sinh  VsTa  L  +  cosh  Vs/a  L 
6Xp(TlS) 


(29) 


The  inverse  transform  must  be  taken.  A  second-order  pole  exists 
at  s  a  0,  and  an  infinite  number  of  simple,  poles  exist  at  the  roots  of 


y/s/a  sinh  Vs/a  L  =  -  cosh  Vs/a  L 


(30) 


or  cotA  =  j-j-A  where  A  =  i y/s/a  L. 

One  can  see  from  Eq.  29  that  the  inverse  transform  will  consist 
of  two  parts,  one  the  inverse  transform  of 


cosh  y/s/q  (L  -  X) 


s(k/h  y/~s/a  sinh  y/sTa"  L  +  cosh  y/sJa  L) 


(31) 


the  other  part  will  have  a  similar  inverse  transform  except  that  it 
will  be  translated  along  the  time  axis  bv  units. 

First,  the  inverse  transform  of 


cosh  Vs /a  (L  ~  X) 


s(k/h  y/s/a  sinh  y/s/a-  L  +  cosh  V s/a  L) 


(32) 
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will  be  located  and  then  integrated  from  0  <  T  <.  T. 


The  inverse  transform  may  be  calculated  by  the  method  of  residues, 
residues  of  the  simple  poles  mav  be  derived  from 

„  .  v*  !W  _/.  „\ 


y  -pk  **?(*  t) 

k  q  (s»)  1  "  1 


where  is  a  ratio  of  polynomials. 


q’  =  s  ~  Js/ct  sinh  \Js/ a  L  +  cosh  sfs/o.  L 

as  h 


+  \Js/a  sinh  y/s/a  L  +  cosh  \fsja  L 


The  second  term  is  zero  since  the  two  parts  of  the  second  term  are 
equated  to  determine  the  roots. 

Performing  the  differentiation  indicated  in  the  first  term  of  Eq.  34 
yields 


[-(l)2  t +  (sc +  { 


sinh  s/s /a  L 


and  changing  to  trignometric  functions  gives 


Substitution  of 


A  =  i  \JsJa  L 

into  the  preceeding  expression  yields 


- 1  (hi.) 


*2  +  fe  +  1  sl"x 


The  numerator  of  Eq.  33,  P,  mav  be  changed  to  trignometric  functions 


— —  a-—-  ^ 

P  cosh  yj s/a  X  =  cos  i  ,/s/a  X  =  cos  —  . 
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Substitution  of  the  above  expressions  into  Eq.  33  yields 


cos  A  X/L  exp(-A  2  aT/L2) 


sin  A 


which  may  be  simplified  to 


Ip  =  -  4 


00 

£ 


cos  A  X/L  sin  X  exp  f-X  ^  aT/L^  ) 
n  n  r  \  n  / 


2A  +  sin  2  X 
n  n 


The  right  side  of  Eq.  36  must  be  integrated  between  the  limits 
of  zero  and  T.  Assume  the  order  of  integration  and  summation  may  be 
interchanged. 

Esin  X  cos  X  X/L  /*T  /  o  9\ 

rLo  expK  d6/L)d6 


00 

£ 


sin  X  cos  X  X/L 
n  n 


„  X/L  /  -  ?v 

- - r  exp(-X  *  aT/LZj  -1  (37) 

n  2  X  \  \  n  / 


*  “T*  A  2(2  1  +  sin  2  X  j  r\  "  ) 

n=l  n  \  n  n J 

The  residue  of  the  second-order  pole  at  s  =  0  must  be  determined. 


P  =  s-K)  I?  "  °^2  T  ^X,s^  exp 


3  [ cosh  \/s/a  (L-X)  exp(sT) 


3s  I  k 


v/s/a  sinh  y/s/a  L  +  cosh  y/sja  L 


=  y/s/a  sinh  y/s/a  L  +  cosh  yjs/a  L^ 

•  sinh  yfsja  (L-X)  exp(sT)  +  T  exp(sT) 

cosh  yj s/a  (L-X)  -  cosh -/sja  (L-X)  exp(sT) 
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f 


1 


f 

[' 


r 


( 


~  ~  cosh  JsJa  L  +  ~  ■_  1  sinh  JsJa  L  +  zr~k= 

h  v  h  2  as  2yas 


2a  h 
•  sinh 


x/s7a  L  j 

all  divided  by|~  yjs/a  sinh  ^ /s/a  L  +  cosh  ^/sTa  hj2  . 

The  limit  of  Eq.  28  as  s  ->0  gives  p  =  T 

The  summation  of  all  residues  yields  the  inverse  transform: 

fexp^-x^2  qT/h2^-l|  sinh  XR  cosh  XR  X/L 


Cf  f  +  4L2 
t  =  —  <  T  + - 


T- 


a 


E 

n=l 


^2X  +  sin  2  X  \ 

\  n  n; 


(38) 


(39) 
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